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ABSTRACT 

Analysis  of  digital  remote  sensing  images  can,  and  should, 
employ  all  the  information  that  is  available  about  the  area 
of  interest.  Maps  often  do  exist,  s.o  that  the  question  arises 
how  map  information  can  be  automatically  used  do  direct  image 
analysis.  With  an  appropriate  organisation  of  digital  map  data 
banks  it  becomes  possible  to  merge  image  and  map.  So  far  we 
have  concentrated  on  areal  features .  We  describe  in  the  report 
the  algorithms,  data  structures  used  for  the  task  and  we  report 
on  our  experiences  using  specific  examples  of  Landsat  images  in 
Southern  Germany  and  Austria  to  automatically  recognize  features 
for  subsequent  rectification. 


RESUME 


L' analyse  des  images  numeriques  de  telfedetection  doit 
utiliser  toute  1 ' inf ormation  qui  existe  d'une  region.  Les 
cartes  g£ographiques  normalement  existent;  il  se  pose  la 
question  comment  utiliser  ces  cartes,  d'une  mani£re  automatique, 
pour  diriger  1' analyse  d'une image.  Une  banque  de  donn6es  numeri¬ 
ques  doit  etre  arrang^e  d'une  manidre  specif ique  pour  relier 
image  et  carte.  Jusqu'd  pr&sent  tous  nos  efforts  se  concen- 
traient  aux  objets  ardaux.  Dans  le  rapport  nous  decrivons  algo- 
rithmes,  la  structure  des  donndes,  et  nous  discutons  les  ex¬ 
periences  obtenues  en  reconnaissant  automatiquement  des  objets 
dans  des scenes  Landsat,  situde  au  sud  de  l'Allemagne  et  en 
Autriche,  et  son  application  au  redressement  des  images. 


ZUSAMMENFASSONG 


Die  Auswertung  digitaler  Fernerkundungsbilder  soil  auf 
der  gesamten  Information  beruhen,  welche  viber  das  Interessensge- 
biet  zur  Verfugung  steht.  Meist  bestehen  Karten,  sodafi  sich  die 
Frage  erhebt,  wie  die  Kartendaten  fUr  eine  Steuerung  der  Bild- 
analyse  in  automatischer  Weise  verwendet  werden  konnen.  Mit  ei- 
ner  geeigneten  Organisationsform  ftir  die  digitalisierten  Karten- 
inhaite  in  einer  Kartendatenbank  wird  es  moglich,  Karte  und  Bild 
zu  verschmelzen.  Bisher  haben  wir  uns  auf  flSchenhafte  Merkmale 
in  Bild  und  Karte  konzentriert,  urn  entsprechend  Bild  und  Karten- 
merkmale  automatisch  aufzufinden.  Wir  beschreiben  in  der  vorlie- 
genden  Arbeit  unsere  Methoden,  Datenstrukturen  und  Erfahrungen 
mit  Beispielen  von  Landsat-Aufnahmen  in  Siid-Deutschland  und  Oster- 
reich  in  der  automatischen  Erkennung  von  Merkmalen  zur  nachfol- 
genden  Entzerrung. 


1,  INTRODUCTION 


The  analysis  of  remote  sensing  imagery  should  rely  on 
available  data  from  the  investigated  area. Generally  this  will 
be  so  called  ground  truth  in  the  form  of  observations  made  in 
the  field  or  of  imagery  of  higher  diagnostic  power,  e.g.  of  a 
larger  scale.  In  many  instances,  however,  remote  sensing  ima¬ 
gery  is  being  analyzed  of  areas  where  excellent  maps  already 
exist.  A  number  of  investigators  have  thus  come  up  with  proce¬ 
dures  to  merge  digital  imagery  with  rhe  existing  digitized  map 
so  that  the  interpreter/user  can  take  full  advantage  of  both 
data  sources.  So  far,  these  procedures  were  essentially  manual, 
with  the  purpose  of  merely  presenting  a  combined  data  set 
(BRYANT  and  ZOBRIST,  1977) . 

The  question  arises  naturally,  how  a  procedure  must  oper¬ 
ate  that  merges  digital  images  and  maps  automatically  and  how 
the  data  present  in  the  map  can  help  to  automatically  analyze 
a  digital  image. 

We  have  taken  it  upon  us  to  study  these  questions.  Obvious¬ 
ly  they  address  a  very  complex  area  that  falls  into  the  catego¬ 
ry  of  automatic  photo-interpretation.  We  are  well  aware  of  the 
fairly  pessimistic  opinions  sometimes  expressed  on  the  prospects 
of  automation  of  photo-interpretation.  Nonetheless  do  we  believe 
that  it  makes  sense  to  support  human  analysis  of  digital  ima¬ 
ges  by  automation  and  to  incorporate  an  available  map  into  this 
process.  There  are  certain  monitoring  tasks  such  as  for  example 
in  the  detection  of  changes  that  could  be  nearly  fully  automa¬ 
ted  in  this  way;  for  other  remote  sensing  applications  this 
may  not  be  the  case;  the  benefits  of  automation  for  those  appli¬ 
cations  may  consist  of  greater  efficiency  of  the  interaction 
between  man  and  machine. 

Map-guidance  in  the  image  analysis  is  central  in  our  effort. 
The  term  map-guided  automatic  photo-interpreation  has  previous¬ 
ly  been  used  by  BARROW  et  al.  (1977).  Essentially  we  see  in  our 
work  three  areas  of  emphasis; 

-  the  organisation  of  the  digital  map  data; 

-  the  analysis  of  areal  features ; 

-  the  analysis  of  linear  features. 

The  tools  required  encompass  a  range  of  manipulating  digital 
graphics  data,  of  image  processing  and  of  automatic  pattern  re¬ 
cognition.  Our  general  approach  to  automatically  merge  image  and 
map  was  described  in  a  previous  paper  (KROPATSCH  and  LEBERL, 

1978  a) .  We  will  review  this  concept  in  the  next  section.  This 
will  be  followed  by  a  description  of  the  manipulation  of  map 
data  (chapters  3  and  4) ,  of  a  discussion  of  image  pre-processing 
(chapter  5),  and  the  establishment  of  relationships  between  maps 
and  image  (chapter  6) . 

In  our  study  we  employ  the  techniques  of  automated  feature 
recognition  to  the  geometric  image  rectification.  A  description 
of  these  fairly  common  procedures  follows.  Finally,  we  go  into 
our  results  as  obtained  with  Landsat  digital  satellite  images . 
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Previous  work  in  this  area  has  been  performed  by  STIES 
et  al.  (1977) ;  they  described  a  comprehensive  image  inter¬ 
pretation  system  that  includes  a  map  data  bank  suitable  for 
map-guided  image  analysis.  BARROW  et  al.  (1977)  and  TENENBAUM 
et  al.  (1978)  have  mainly  concentrated  on  the  analysis  of 
linear  features  in  digital  images.  With  satellite  images  at 
small  scales,  however,  it  is  promising  to  study  mainly  areal 
features.  After  all,  satellite  remote  sensing  draws  its  essen¬ 
tial  motivation  from  the  analysis  of  natural  surface  features: 
it  is  for  this  reason  that  scales  can  remain  comparatively 
small,  while  the  study  of  manmade  linear  objects  would  require 
larger  scales  and  thus  airborne  remote  sensing. 

However,  some  results  have  been  obtained  also  with  linear 
features  from  satellite  scanner  data. 

It  should  be  emphasized,  that  the  techniques  may  have  a 
broad  spectrum  of  uses ,  not  only  for  rectification  but  also 
for  other  tasks.  As  an  example  one  may  point  to  the  work  of 
FLONZAT  et  al.  (1979)  who  based  multi-spectral  classification 
on  training  areas  which  are  taken  from  a  map  and  transformed 
into  the  image. 

If  we  go  back  somewhat  we  find  that  an  early  automatic 
image  recognition  application  was  the  measurement  of  pricked 
points  (KREILING,  1976)  and  of  reseau  points  (ROOS,  1975)  . 

The  associative,  complex  mental  interpretation  of  an  ima¬ 
ge  by  the  human  expert  will  hardly  ever  he  matched  by  the 
computer.  However,  specific  tasks  may  be  ideally  suited  for 
automation.  In  the  current  context  we  are  certainly  not  tack¬ 
ling  the  most  trivial  application  but  one  which  in  its  com¬ 
plexity  helps  one  to  understand  the  complexity  of  the  problem 
in  general.  We  may  take  the  statement  of  a  prominent  pattern 
recognizer  for  whatever  it  is  worth,  but  ALEKSANDER  (1978) 
claimed, that  "most  problems  of  automated  pattern  recognition 
have  solutions;  difficulties  only  exist  in  the  speed  and  ge¬ 
neral  applicability  of  these  solutions".  Our  current  study  is 
aimed  at  the  understanding  of  solutions  to  the  automated  re¬ 
cognition  of  map  features  in  images.  Our  first  application  is 
to  use  these  features  as  control  points  for  a  geometric  recti¬ 
fication. 
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2.  GENERAL  CONCEPT  OF  ARSIM 


The  task  of  ARSIM  (ARSIM  is  an  acronym  for  Automatic 
Registration  of  Satellite  Images  and  Maps)consits  of  re¬ 
gistering  a  given  satellite  image  with  the  corresponding 
map.  Because  of  local  distortions  of  the  image  the  re¬ 
gistration  cannot  be  done  with  a  simple  global  transforma¬ 
tion  defined  a  priori.  A  common  way  to  register  images 
with  each  other  is  the  use  of  control  points .  This  method 
is  also  suitable  for  the  present  task.  It  has  first  to 
produce  enough  control  points  so  that  one  can  rectify  the 
image  geometrically.  Figure  2.1  presents  a  general  program- 
flow  and  shows  that 


(1) 

program  ARSIM 

(2) 

begin 

(3) 

REGISTER  (map,  image,  control  points) ; 

(4) 

if  sufficient  control  points  then 

(5) 

RECTIFY  (image,  control  points. 

new  image) 

(6) 

else 

(7) 

print  ( 1  map  and  image  cannot  be 

registered' ) ; 

(8) 

end. 

Figure  2.1:  Main  program  for  project  ARSIM. 


the  problem  splits  into  two  mainly  independent  procedures 
REGISTER,  which  registers  map  and  image  delivering  control 
points,  and  RECTIFY,  which  rectifies  the  geometry  of  the 
image  with  the  aid  of  the  control  points  to  result  in  a  new 
image.  The  new  image  has  then  the  geometry  of  the  map.  The 
method  of  operation  is  described  in  detail  in  chapter  7. 

Before  the  processing  of  the  lower  level  procedures  is 
described  in  more  detail  an  overview  of  the  registration 
process  is  discussed  by  means  of  a  flow-diagram  (figure  2.2) . 

In  this  flowchart  processing  steps  (1)  to  (9)  correspond 
to  the  procedure  REGISTER,  while  step  (10)  correspond  to 
procedure  RECTIFY.  The  following  chapters  3  through  7  describe 
the  functions  of  the  lower  level  procedures  in  detail  mostly 
using  a  general  purpose  Problem  Description  Language  (PDL) . 


2.1  THE  REGISTRATION  PROCESS 


Following  NACK  (1975)  we  define  in  the  digital  image  a 
square  (rectangular)  Area  of  Interest  (AI)  that  is  being 
processed  at  any  given  time. 


i 


For  processing-time  considerations  the  digital  image  to 
be  registered  with  the  map  devides  into  a  number  of  AIs.  These 
are  easier  to  handle  than  the  total  image.  They  restrict  the 
search  and  speed  up  the  procedure.  Thiu  logical  processing 
unit  subdivides  further  the  operations  of  procedure  REGISTER 
into 

(A)  operations  with  AIs  (steps  (3)  until  (6)  in  figure  2.2), 
and 

(B)  operations  within  AIs  (steps  (7)  until  (9)  in  fig.  2.2). 

Operations  (A)  form  procedure  REGISTER  (figure  2.3),  the 
operations  from  (B)  are  collected  in  a  subroutine  called 
OVERLAY  and  described  in  chapter  2.2. 


(1)  procedure  REGISTER  (map,  image,  controlpoints) ; 

(2)  begin 

(3)  controlpoints :=  .empty; 

(4)  repeat 

(5)  begin 

(6)  SELECT  (AI) ;  0  input  an  Area  of  Interest  0; 

(7)  TRANSFORMWINDOW  (AI,  MW) ; 

(8)  MAPLOAD  (map,  MW) ; 

(9)  PREPROCESS  (image,  AI) ; 

(10)  OVERLAY  (image,  AI,  map,  MW,  controlpoints  AI) ; 

(11)  controlpoints:  =  controlpoints+controlpoints  AI  ; 

(12)  end  until  no  input  left; 

(13)  end . 


Figure  2.3:  Operations  with  Areas  of  Interest  (AI) . 


In  step  6  (figure  2.3)  we  start  from  a  digital  map  base  that 
is  available  in  a  format  compatible  with  the  image  processing 
system  DIBAG,  which  contains  the  satellite  image  to  be  re¬ 
gistered.  SELECT  (AI)  reads  an  AI,  which  can  be  chosen  manu¬ 
ally  by  the  following  criteria: 

-  it  has  a  predefined  window  width, 

-  in  that  window  there  is  a  certain  amount  of  patterns, 
which  can  be  optically  distinguished  in  the  image, 

-  for  most  of  these  patters  in  the  image  exist  corres¬ 
ponding  objects  in  the  map. 

Each  pixel  of  a  digital  satellite  image  can  be  located  in  a 
map  projection  with  a  certain  positional  accuracy.  For  LANDSAT, 
maximum  errors  have  been  demonstrated  to  amount  to  not  more 
than  5  km  (COLVOCORESSES  et  al . ,  1973).  Data  used  to  achieve 
this  accuracy  are  the  satellite  position  and  attitude  values 
delivered  with  each  image.  More  recently,  orbit  parameters  in 
LANDSAT  permit  an  even  better  accuracy  to  be  achieved  without 
ground  control  data,  namely  about  10  pixel  diameters  or  800  m. 
Using  the  map  locations  of  the  corner  points  of  AI  we  define 
the  corresponding  map  window  (MW) (procedure  TRANSFORMWINDOW 
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in  step  (7)  of  figure  2.3)  .  In  order  to  contain  all  of  AI, 
the  window  has  to  be  larger  than  AI  due  to  the  limited  po¬ 
sitional  accuracy  (compare  figure  2.4). 

Then  MAPLOAD  loads  a  portion  of  the  map  data  bank  which 
comprises  only  objects  completely  in  MW  (step  8  of  figure  2.3) . 
The  structure  of  the  data  bank  is  described  in  chapter  3 . 

Depending  upon  the  types  of  objects  contained  in  that 
portion  of  the  map  we  are  now  able  to  call  routine  PREPROCESS 
(step  9  of  figure  2.3)  for  the  AI  in  the  image.  For  each  ob¬ 
ject-type  there  exists  a  special  method  of  pre-processing  of 
the  image  data  (for  more  detail  see  chapter  5), so  that  when 
we  search  for  some  object  we  do  it  on  an  image  that  is  pre- 
processed  for  this  object-type. 


2.2  REGISTRATION  IN  AN  AREA  OF  INTEREST. 


After  all  the  preparations  and  restrictions  we  are  now 
able  to  turn  to  the  main  task:  identification  of  objects 
from  the  map  in  the  image.  This  is  done  in  procedure  OVERLAY 
and  consists  of  5  major  operations  as  shown  in  figure  2.5: 


-  the  selection  of  an  object  from  MW  ( OB JECTSELECTION ) ; 

-  the  extraction  of  the  object  from  the  data  bank  and 
its  transformation  in  a  format  which  allows  a  compari¬ 
son  with  the  image  (MAPTRANSFORM) ? 

-  the  recognition  of  the  object  in  the  image  (RECOGNIZE) ; 

-  the  VERIFICATION  that  the  pattern  found  in  the  image 
corresponds  to  the  original  object  and 

-  the  search  for  controlpoints  (CONTROLP) . 


(1)  procedure  OVERLAY  (image,  AI,  map,  MW,  controlpoints) ; 

(2)  begin 

(3)  controlpoints:  =  empty; 

(4)  OB JECTSELECTION  (map,  MW,  list); 

(5)  while  list  not  empty  do 

(6)  begin 

(7)  object:  =  first  (list); 

(8)  list:  =  list  -  first  (list) ; 

(9)  MAPTRANSFORM  (map,  object,  binarymatrix) ; 

(10)  RECOGNIZE  (image,  AI,  binarymatrix,  pattern); 

(11)  if  VERIFICATION  (object,  image,  pattern)>accepttresh  then 

(12)  CONTROLP  (object,  pattern,  controlpoints) 

(13)  end; 

(14)  end. 


Figure  2.5:  Overlay  of  image  AI  and  MW. 


Features  within  MW  are  sorted  according  to  their  probabi¬ 
lity  of  being  identifiable  in  the  satellite  image.  This  probabi¬ 
lity  can  be  (automatically)  adjusted  as  the  result  of  previous 


experiences.  Consequently,  the  procedure  has  a  learning  ca¬ 
pability.  The  procedure  OBJECTSELECTION  (figure  2.6)  provi¬ 
des  a  sorted  list  of  these  features.  The  most  probable  fea¬ 
ture  within  MW  is  selected  and  a  search  begins  in  the  image-AI 
to  find  the  corresponding  detail.  MAPTRANSFORM  (see  figure 
4.2)  extracts  a  boundary  polygon  of  this  "object"  from  the 
map  data  bank,  it  linearly  transforms  the  polygon  into  the 
image  coordinates,  using  the  same  function  as  in  TRANSFORM- 
WINDOW  and  builds  a  binary  matrix  which  determines  if  a  pixel 
of  the  AI  is  either  inside  the  polygon  or  outside.  The  binary 
matrix  of  the  object  serves  as  input  to  procedure  RECOGNIZE 
(figure  2.7).  It  breaks  down  the  recognition  process  into 
4  different  methods,  which  are  described  in  chapter  6.  The 
result  of  the  search  is  another  binary  matrix  called  "pat¬ 
tern  " . 


(1)  procedure  OBJECTSELECTION  (map,  MW,  list) ? 

(2)  begin 

(3)  list:  =  empty? 

(4)  for  all  objects  €  MW  do 

(5)  begin 

(6)  recogn:  =  probability  (object) ? 

(7)  if  recogn  > recogthresh  then  list:=  list  + (object , recogn) 

(8)  end? 

(9)  SORT  (list) 

(lO)  end. 


Figure  2.6:  Selection  and  Sort  of  objects  in  MW. 


Next  the  identity  of  the  map  and  image  feature  has  to  be 
verified  by  a  comparison  of  the  "object"  and  the  "pattern" 
(VERIFICATION) .A  verification  measure  is  compared  to  a  certain 
"accepttresh"  to  include  the  worst  case:  a  reject. 

A  last  step  (12  in  figure  2.5)  is  the  search  for  singul¬ 
arities  of  the  feature  so  that  unique  homologue  pairs  of  co¬ 
ordinates  can  be  identified  in  both  the  image  and  the  map 
(CONTROLP) .  As  this  is  an  essential  concept  for  the  rectifi¬ 
cation  process,  chapter  7  contains  a  description  of  it. 


(1)  procedure  RECOGNIZE  (image,  AI,  binarymatrix,  pattern) 

(2)  begin 

(3)  case  search  method  (object)  of 

(4)  begin 

(5)  1:  SHIFT  (image  ,  binarymatrix,  pattern) 

(6)  2:  THRESH  (image,  binarymatrix,  pattern) 

(7)  3:  ADAPT  (image,  binarymatrix,  pattern,  D) 

(8)  4:  LINDET  (image,  binarymatrix,  start  (object), 

goal  (object) ,  pattern) ? 

(9)  end 

(10)  end. 


Figure  2.7:  Recognizing  an  object  by  different  methods. 
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3.  OPERATIONS  WITH  MAP  DATA 


3.1  GENERAL 


Processing  of  digital  map  data  is  basically  different 
from  that  of  images;  we  are  not  dealing  here  with  raster 
type  data,  but  with  vectors.  The  area  of  work  is  often  times 
denoted  by  the  expression  "computer  graphics"  and  is  seen 
to  be  the  fastest  growing  portion  of  the  computer  industry. 

In  the  current  context  we  are  dealing  with  computer 
graphics  in  a  specific  way:  we  are  processing  vector  data 
and  raster  information  and  need  conversion  of  one  into  the 
other.  This  is  not  a  trivial  task,  even  if  it  may  seem  so 
at  superficial  inspection. 

Graphical  data  banks  consist  of  data  and  a  computer 
program  system  for  manipulation  of  the  data.  A  vast  body  of 
literature  exists  in  this  area,  also  from  the  particular 
point  of  view  of  cartography  (e.g.  FRANK,  1979.;  WEBER,  1978). 
Depending  on  specific  requirements  and  uses  data  banks  may 
be  simple  or  very  complex.  Data  structures  define  complexity 
and  in  turn  result  from  the  type  of  relationships  that  must 
be  defined  among  elements  of  the  data  bank. 

Data  structures  can  roughly  be  grouped  into  two  groups: 

(a)  one  without  explicit  definition  of  relationship  among 
data  elements  (e.g.  object  oriented  sequential  or 
raster  data  structure)  ,* 

(b)  one  with  explicit  definition  of  relationship  (topology, 
neighbourhood,  positional  reference) . 

Most  of  currently  implemented  cartographic  data  banks 
belong  to  category  (a)  without  a  thorough  definition  of  re¬ 
lationship  among  data  elements . 

In  our  current  application  we  face  the  task  of  defining 
a  map  data  bank  of  category  (b) .  It  must  thus  be  possible 
to  establish  very  quickly  interconnections,  neighbourhoods, 
relationship  among  objects  of  a  map,  much  like  the  process 
of  a  visual  inspection  of  map  contents  by  a  human  interpreter. 
A  topological  data  structure  will  thus  be  required.  It  is 
straight-forward  to  define  a  set  of  functions  to  be  satis¬ 
fied  by  the  data  bank: 

-  selection  and  extraction  of  single  elements  (objects) 
of  the  data  bank,  using  names,  properties,  positions; 

-  establishment  of  connections  between  individual  ob¬ 
jects  in  geometrical  and  thematic  set  operations; 

-  computations  with  elements ;  . 

-  generation  of  lists  and  printouts; 

-  graphical  output  in  vector  or  raster  format. 

This  enumeration  may  not  be  complete  but  represents  a  set 
of  important  features  to  be  had  by  the  data  bank  in  its 
application  to  work  with  maps  and  images . 
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We  will  describe  in  this  chapter  the  basic  data  struc¬ 
tures,  algorithms  and  uses  of  the  data  bank.  Detail  is 
appropriate  to  instruct  a  prospective  user  of  the  data  bank. 
However,  it  needs  to  be  stressed  that  the  data  bank  is  experi¬ 
mental  in  nature,  and  the  current  implementation  is  motivated 
solely  be  the  need  to  experiment  with  combined  sets  of  ima¬ 
ge  and  map  data.  Numerous  improvements  will  have  to  be  made 
in  the  areas  of  user  friendly  data  acquisition  and  editing. 


3.2  DATA  STRUCTURES 


3.2.1  General 


We  deal  with  data  that  are  available  in  the  form  of 
strings  of  x,y-coordinates  representing  lines  with  a  speci¬ 
ally  identified  beginning  point.  It  is  not  relevant  here 
how  these  points  were  digitized  -  be  it  manual  digitizing, 
line  following  or  raster  scanning  with  subsequent  vectorisa  .- 
tion. 


There  exists  a  multitude  of  data  structures  for  various 
types  of  graphical  data  banks.  Methods  to  store  those  features 
of  a  map  in  digital  form  that  are  relecant  in  our  current 
problem  could  be: 

(i)  Each  feature  is  represented  by  its  contour  polygon 
(object-oriented,  sequential  data  structure) ; 

(iii)  Each  feature  is  defined  by  a  sequence  of  line  segments. 
Each  segment  separates  exactly  2  features  and  consists 
itself  of  a  polygonal  sequence  of  points.  The  map  re¬ 
presents  a  "planar  graph". 

(iii)  Raster  presentation; 

(iv)  Geo-coded  or  positionally  defined  data,  with  the  special 
case  of  the  quad  tree^  (WEBER,  1978;  DYER  et  al.,  1979); 

2) 

(v)  Description  using  the  Freeman-Chain  Code  (FREEMAN, 
1979). 

Since  we  require  a  connection  between  map  data  and  digital 
images  we  need  to  stress  that  these  images  exist  in  the  form 
of  matrices  of  gray  values.  This  corresponds  to  the  concept 
of  a  raster  structure  for  graphical  data  as  used  by  WEBER 
(1978),  BRtfGGEMANN  (1978),  HARRIS  (.1979)  ,  HARRIS  and  PRESTON 
(1979),  GOODCHILD  (1979)  and  others.  It  would  thus  be  logical 
to  consider  structure  (iii)  for  a  data  bank. 


i)  "Quad  Tree"  is  the  name  of  a  tree  structure  of  4th  order. 

Its  root  is  a  square  of  the  base  plane  containing  the  entire 
feature.  There  exist  four  "sons"  of  a  "node"  by  subdividing 
the  square  in  4  equal  segments.  In  each  "leave"  of  the  tree 
(i.e.  in  a  node  without  successor)  one  finds  whether  its 
square  belongs  to  the  basic  square  in  the  base  plane  or  not. 


This  is  based  on  a  raster  representation.  It  describes  the 
contour  of  a  feature  or  object  using  a  sequence  of  small 
vectors  (standardized  and  coded)  which  connect  adjacent 
raster  points. 


2) 
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However,  with  the  map  features  In  a  raster  format  one 
encounters  Immediately  problems  due  to  storage,  redundancy, 
differences  of  raster  mesh  size  and  orientation  and  complexity 
of  required  transformation  routines. 

If  we  evaluate  the  five  classes  of  data  structures  we 
find  that  an  object-oriented,  sequential  method  Ci)  does  not 
compare  well  with  method  Cii)  because  of  redundancy:  all  lines 
must  be  digitized  and  stored  twice. 

Reference  to  position  in  method  Civ)  appeals  due  to  its 
efficiency  when  addressing  points  or  features  in  a  two- 
dimensional  space.  Judging  method  (v)  shows  that  this  is 
suitable  particularly  in  cases  where  data  were  digitized  in 
a  raster  format  by  a  scanner.  Current  assumptions  are,  how¬ 
ever,  that  data  exist  in  the  form  of  vectors.  But  method  (v) , 
the  Freeman-Chain  Code,  still  has  an  advantage  of  reduced 
storage  requirements.  This  is  set  off  by  increased  difficult¬ 
ies  for  geometric  transformations . 

These  various  data  structures  are  not  mutually  exclusive. 
A  feature-oriented,  topologically  meaningful  structure  acc. 
to  (ii)  can  exist  along  with  a  positional  one  such  as  (iv) . 
This  then  would  permit  one  to  get  an  answer  to  two  questions 
like: 

-  Where  is  a  specific  feature  and  what  are  its  neighbours? 

-  What  feature  exists  at  a  specific  position? 


3.2.2  A  structure  to  merge  map  and  image 


Prior  to  a  decision  on  the  data  structure  one  needs  to 
identify  the  principles  of  procedures  to  merge  map  and  image. 
We  do  see  two  main  alternatives: 

(a)  A  synthetic  image  is  generated  from  the  graphical  data; 
actual  and  synthetic  images  are  then  conventionally 
cross-correlated . 

(b)  Map  data  is  structured  into  elementary  information, 
namely  features  or  objets.  Of  these  one  selects  specific 
ones  for  a  meaningful  connection  with  a  given  image . 

The  first  route  must  fail  first  of  all  when  maps  contain 
features  that  do  not  appear  in  an  image:  political  boundaries, 
for  example,  or  in  the  reserve  case,  clouds;  furtheron  if 
map  features  are  generalized  as  with  roads.  We  must  realize 
that  maps  have  primarely  been  produced  to  portray  the  earth 1 s 
surface  for  interpretation  by  a  human.  We  therefore  turn  our 
attention  preferably  to  the  second  method  (b)  above. 

One  not  only  has  available  for  a  merging  process  the 
shape  and  position  of  an  object,  but  also  its  topological 
environment.  A  rock  area,  for  example,  is  more  difficult  to 
identify  in  an  image  if  it  is  surrounded  by  glaciers  as 
compared  to  forests. 
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This  leads  then  to  the  use  of  a  topological  data  structure 
(WEBER,  1978)  as  the  highest  level  of  the  data  bank,  accord¬ 
ing  to  method  (ii)  in  section  3.2.1.  This  structure  will  be 
discussed  in  detail  later.  The  highest  level  of  the  structure 
presents  the  position  of  features  with  respect  to  one  another. 
There  still  remains  the  task  to  store  thematic  information 
and  to  exactly  position  the  feature  (geometric  information) . 
Since  the  topological  relationships  have  been  defined  already 
in  the  highest  level  one  can  make  use  of  the  simplest  data 
structure  for  the  lower  levels;  this  is  the  sequential  one. 

It  is  here  that  features  are  described  geometrically  and 
thematically. 

A  data  bank  is  thus  organized  in  three  levels  of  infor¬ 
mation  (Fig.  3.1): 

-  topological  relationship  (structure)  ; 

-  sequential  lists  of  objects,  i.e.  of  features 
(regions)  and  of  lines  (region  list,  line  list) ; 

-  polygons  stored  as  a  sequential  coordinate  list 
(coordinate  file) . 

The  next  three  sections  will  present  the  data  bank, 
followed  by  two  data  structures  useful  for  the  current 
problem:  the  coordinate  list  of  an  ob j ectfe  contour,  and  the 
binary  matrix.  The  former  serves  as  an  intermediate  storage 
and  for  transformation  of  coordinate  systems;  the  latter  is 
needed  for  geometric  set  operations. 

This  structure  also  permits  search  processes  based  on  a 
positional  reference,  however,  in  an  inefficient  manner.  For 
improved  economy  of  this  type  of  search  processes  one  would 
have  to  define  the  positional  reference  set  over  the  topo¬ 
logical  level. 


3.2.3  The  topological  structure,  the  graph. 


In  a  formal  manner  one  can  interpret  a  map  to  be  a 
planar  graph*) : 


G  *  <R,L) 

R  ...  set  of  regions  (features,  nodes); 

L  ...  set  of  lines  (edges),  separating  2  regions. 


l)  Definition:  a  (finite)  graph  G  =  (X,U)  is  composed  of  a 
pair,  consisting  of  the  finite  set  X,  the  "nodes",  and  a 
set  U  c  X  x  X,  the  "edges".  A  graph  is  "planar"  if  it  is 
possible  to  represent  it  in  a  plane  in  such  a  way  that 
all  nodes  are  separable  and  that  all  edges  are  simple 
curves,  without  intersections  or  double  points,  and  that 
two  edges  only  meet  at  the  end  points  ( SAKAROVITCH ,  1975)  . 
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The  graph  should  be  stored  digitally  on  a  minimum  of 
storage  since  all  inputs/outputs  of  data  are  based  on  this 
structure.  This  can  be  achieved  by  a  reduction  and  compress¬ 
ion  of  the  information  kept  in  the  smallest  storage  units, 
so  called  items .  Information  is  reduced  to  pointers  describing 
topological  relationships.  There  are  3  types  of  items  (REG, 
LIN,  LIST) : 

REG:  =  (  @  pointer  RA  to  a  region  list;  @  pointer 

PLIST  to  a  LIST-item) . 


LIN:  =  (  code  262143 n  ;  @  pointer  LA  to  the  line 

list;  ®  pointer  PR1;  @  pointer  PR2) * 

LIST:  =  (  pointer  PL  to  a  LIN-item;  (2)  pointer 

PNEXT  to  the  next  LIST-item) . 


Pointers  PRl  and  PR2  point  to  two  REG-i terns  which 
correspond  to  the  regions  separated  by  a  line. 

On  a  36-bit-computer  each  REG-item  exists  of  1  computer 
word  (2  halfwords) ,  each  LIN-item  of  2  words  (  4  halfwords) 
and  each  LIST-item  of  1  word  (2  halfwords) . 

The  mode  of  data  acquisition  and  properties  of  features 
lead  to  definition  of  two  types  of  regions: 

-  areas  and 

-  skeletons . 

Skeletons  are  areas  represented  by  lines  such  as  rivers 
or  roads.  Two  reasons  exist  for  this  separate  definition  of 
skeletons : 

(a)  they  permit  digitizing  of  such  linearly  extended  areas 
as  if  they  were  lines;  there  is  no  need  to  digitize  all 
contours  of  such  an  area  -  this  would  be  a  redundant 
effort; 

(b)  the  width  of  skeletons  can  be  selected  at  a  later  stage  - 
a  feature  of  importance  when  various  scales  are  needed 

in  visualizing  the  map  contents. 

Skeletons  differ  from  areas  in  the  graph:  the  pointer 
LA  of  a  border  line  is  replaced  by  a  specific  code  (the  num¬ 
ber  262142  in  our  case) .  This  can  be  done  because  in  the  case 
of  a  skeleton  there  are  no  specifically  digitized  border  lines. 


1 )  Redundant  code  to  fill  up  2  computer  words. 
Denotes  a  LIN-item. 


Figure 


L 


3.2:  Digital  map  of  Graz 


data  bank  of  Graz  around  skeleton  S7 
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However,  the  graph  needs  a  LIN-item:  this  is  then  a  pseudo¬ 
line  which  was  not  digitized.  A  skeleton  is  thus  an  area 
with  pseudo-lines  as  contour.  Beginning  and  end-points  of 
a  skeleton  are  two  specific  LIST-items  where  PL-pointers 
are  replaced  by  a  code  for  the  beginning  (the  number  262140 
in  our  case)  and  end  (the  number  26241).  Figure  3.2  illustra¬ 
tes  a  portion  of  a  data  base  in  the  area  of  Graz,  digitized 
off  maps  at  scale  1:50  000.  Elements  of  this  data  base  will 
be  used  to  illustrate  the  topological  data  structure. 


3.2.4  Illustration  of  the  graph  using  a  skeleton. 


Taking  a  specific  skeleton  from  fig.  3.2  we  describe 
its  pertinent  structure  in  figures  3.3  (a) and  3.3(b). 

Fig.  3.3(a)  presents  a  skeleton  extending  between  points 
A  and  E,  and  with  the  name  S7.  As  a  region  or  area  it  is 
delimited  by  a  series  of  pseudo-lines  P313,  P354,  P367, 

...,  P1079.  A  printout  of  the  graph  is  shown  in  fig.  3.3(b), 
with  only  the  data  repeating  to  skeleton  S7 .  Rows  or  lines 
in  the  computer  printout  are  shifted  left-right  depending 
on  whether  it  concerns  a  LIN-,  LIST-  or  REG-item.  Each  line 
of  the  printout  contains  a  line  number,  increasing  regularly 
from  1  to  n,  and  the  content  of  a  36-bit-computer  word  with 
two  integers.  Each  LIN-item  consists  of  2  rows,  each  REG- 
and  LIST-item  of  one  row. 

To  explain  the  printout  in  more  detail.,  let  us  take 
the  first  item  belonging  to  skeleton  S7.  It  is  in  rows  313 
and  314  of  the  graph,  with  the  following  format: 

Row  number  313:  T262143  LA  “1 

Row  number  314:  LpRl  PR2J 


This  is  a  LIN-item  as  is  obvious  from  code  262143.  The  value 
of  LA  is  262142;  we  deal  this  with  a  pseudo  line  separating 
two  regions .  The  first  region  has  its  REG-item  in  row  309  as 
decumented  by  the  value  309  of  pointer  PRl,  the  second  re¬ 
gion  in  row  1070  (sea  pointer  PR2) .  Figure  3.3(b)  presents 
all  LIN-items  relating  to  the  second  region  with  REG-item 
in  row  1070.  Other  rows  of  the  printout  were  suppressed.  REG- 
item  1070  has  an  RA-pointer  7  and  therefore  represents  ske¬ 
leton  s7  (compare  the  region  file,  section  3.2.6). 

Then  LIN-items  in  rows 342,  348,  354  separate  the  region 
in  row  1070  from  a  region  in  row  328.  In  the  region  file 
that  region  has  the  name  EGGENBERG. 

LIST-items  serve  to  accelerate  the  search  for  those  LIN- 
items  that  belong  to  a  given  region.  From  row  1070  we  find, 
using  pointer  PLIST,  that  the  first  LIST-item  of  the  region 
EGGENBERG  is  in  row  1071.  Each  LIST-item  consists  of: 
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Figure  3.4(b):  Sections  of  the  graph  around  region  197, 
Eggenberg . 


-25- 


Row  number:  PL  PNEXT. 

Each  pointer  PL  points  to  the  next  LIN-item  or  designates  the 
start  or  end  of  the  contour  of  a  region.  PNEXT  points  to  the 
row  with  the  next  LIST-item. 


3.2.5  Illustration  of  the  graph  using  an  area. 


Figures  3.4(a),  (b)  show  two  areas,  one  enclosed  by  the 
other.  The  pseudo-line  has  special  significance  in  this  case 
that  is  typified  by  an  island  in  a  lake:  non-connected  parts 
of  the  boundary-polygon  (contour)  exist  in  the  form  of  the 
edges  of  the  lake  and  island.  The  pseudo-line  connects  both 
and  must  be  passed  twice. 

From  this  and  the  illustration  of  the  skeleton  in  fig. 3. 3 
it  is  now  straight  forward  to  follow  the  section  of  the  graph 
for  the  areas  EGGENBERG  and  SCHLOSSEGGEN . 

Row  328  contains  the  REG-item  of  region  EGGENBERG  with 
pointer  PLIST  pointing  to  row  331  with  the  first  LIST-item. 
This  in  turn  has  pointer  329  to  the  (first)  LIN-item  in  row 
329.  Among  the  LIST-items  there  are  two  rows  with  pointer  PL 
pointing  to  the  same  LIN-item  in  row  332:  this  designates  a 
pseudo-line  (LA  equals  262142)  separating  area  328  from  328 
or  more  appropriately  passing  through  to  area  328. 

Region  EGGENBERG  is  delimited  by  region  SCHLOSSEGGEN,  by 
skeleton  S44  and  by  the  undefined  background  (S44  internally 
has  the  name  S00000000045)  .  This  is  the  reason  why  figure 
3.4(b)  also  contains  the  graph  sections  for  these  regions. 


3.2.6  Line-  and  region  files 


Line-  and  region  files  form  the  connection  between  the 
graph  and  the  lowest  level  consisting  of  the  sequence  of  points 
in  the  coordinate  list.  Line-  and  region  files  contain  descrip¬ 
tive,  thematic  and  statistical  information  concerning  a  line 
or  region.  An  element  of  these  files  is  composed  as  follows: 

(a)  Name  .  12  characters.  During  data  acquisition  or 

tagging  each  line  and  skeleton  automatical¬ 
ly  are  assigned  names .  For  lines  this  is 
LOOOOOOOOXXX  /  for  skeletons 

SOOOOOOOOXXX  ,  where  'XXX'  represents 

successive  running  numbers .  Names  of 
regions  can  be  freely  selected  by  the 
user  in  the  process  of  building  up  the  graph. 

4  values .  The  first  two  represent  the  small¬ 
est  x,y  coordinates  of  a  line  or  skeleton, 
the  last  two  the  dimension  in  x-  and  y-direct- 
ion  of  the  data  bank  coordinate  system. 


(b)  Windows 


-26- 


(c)  Pointer  . points  to  the  respective  item  in  the 

graph . 

(d)  Coordinates  .  2  values  delimiting  the  range  of  rows 

in  the  coordinate  list  for  the  polygon 
of  a  line  or  skeleton. 

(e)  Starting  point  ...  of  a  line  or  skeleton. 

(f)  End  point  .  of  a  line  or  skeleton.  These  points  per¬ 


mit  rapid  compilation  of  an  areals  bounda¬ 
ry  or  contour  without  directly  address¬ 
ing  the  coordinate  list.  Areas  do  not 
have  start-  and  end  points,  these  exist 
only  for  lines  or  skeletons.  With  areas 
the  values  for  start-  and  end  points 
are  filled  by  zeroes . 

(g)  15  more  positions  .  for  statistical  measures,  for  example 

the  type  of  a  region.  In  our  specific 
application  we  define: 

0  undefined 

1  water  (not  skeleton!) 

2  river,  canal  (only  skeleton!) 

3  built-up  area 

4  range  land  ,  agriculture 

5  forest 

6  rock 

7  road  (skeleton) 

8  rail 

Figure  3.5  presents  an  example  of  a  printout  of  3  parts  of 
a  line-  and  region  file. 


3.2.7  Coordinate  list 


This  file  contains  point  coordinates  forming  polygonal 
lines.  Each  point  can  be  addressed  by  his  position  in  the 
sequential  file.  Figure  3.6  shows  an  example  of  a  printout 
of  the  coordinate  list  for  lines  L00000000031,  LOOOOOOOOQ32, 
LOOOOOOOOO 33. 


3.1* .8  Coordinate  list  of  the  boundary  of  an  area 


We  use  an  intermediary  data  structure  before  and  after 
geometrical  transformation  of  polygons .  Polygons  here  are 
closed  boundaries  of  areas.  The  structure  is  sequential.  Two 
areas  where  one  is  totally  enclosed  by  the  other  are  represented 
as  follows: 

(i )  Each  line  segment  is  preceded  by  a  pair  of  values  marking 
the  start-  and  ending  indices  of  the  sequence  of  points  of 
the  line  segment. 
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Figure  3.6:  Example  of  a  printout  for  coordinate  list  of 
map  Graz  (compare  figure  3.5). 
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:  Coordinate  list  of  a  polygon  describing  area 

Eggenberg.  (*)  marks  value  pairs  preceding  a 
section  of  a  polygon. 
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(ii)  These  line  segments  are  written  sequentially  into  the 
coordinate  list. 

Figure  3.7  illustrates  such  a  coordinate  list  with  the  pre¬ 
ceding  pairs  of  values  here  as  indices  1  and  49. 


3.2.9  Binary  matrix  (raster) 


For  logical  operations  with  graphical  data,  such  as 
operations  of  intersection,  union  and  taking  the  negative 
the  raster  format  is  well  suited.  Therefore  a  provision 
exists  in  the  data  b.ank  to  convert  areas  for  which  a  bound¬ 
ary  polygon  exists  to  a  raster  format.  The  appropriate  data 
structure  is  a  binary  matrix.  Coordinates  are  defined  by  a 
window  preceding  the  matrix.  The  window  definition  is  the 
same  as  that  in  the  line-  and  region  file.  It  is  followed 
by  a  2-dimensional  matrix  M: 

M  (i,j)  =1  if  (i,  j )  £  area 

M  (i,J)  =  0  elsewhere. 


The  elements  of  the  matrix  are  ordered  row  by  row  in  a 
compact  format:  the  last  bit  of  one  row  is  directly  followed 
by  the  first  bit  of  the  next  row.  This  is  currently  implemented 
to  form  arbitrary  rectangular  binary  matrices  with  72  OOO 
elements . 


Figure  3.8  shows  region  EGGENBERG  as  a  binary  matrix. 


3.3  ALGORITHMS  FOR  DATA  INPUT  AND  OUTPUT 


The  current  implementation  of  computer  programs  is  to 
serve  in  the  main  task  of  merging  map  and  image  data.  Limita  - 
tions  of  resources  has  therefore  mainly  affected  the  user  er- 
gonomy  in  data  acquisition  and  editing.  To  further  illustrate 
the  current  implementation  stage  we  outline  the  data  input  and 
output.  Actual  applications  software  now  exists  for  merging 
maps  and  images.  This  is  described  in  a  separate  chapter  4. 


3.3.1  Data  Acquisition 


(a)  Step  1:  Digitizing  of  points,  lines  and  skeletons. 

We  currently  do  off-line  digitization  on  a  manual 
digitizing  table1)  and  then  transfer  the  data  to  a  disk. 


1 )  This  is  possible  due  to  the  cooperation  of  the  Agrar- 
technische  Abteilung  of  the  Government  of  Styria. 


a 
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Figure  3.8:  Printout  of  binary  matrix  of  area  EGGENBERG 
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Editing  is  with  a  text-editor  to  preparate  the  data  for 
establishment  of  the  coordinate  list,  line  file  and  the  ske¬ 
leton  part  of  the  region  file.  This  is  done  from  the  raw  di¬ 
gitizer  data  in  program  MAPREAD.  This  provides  the  following 
elements  in  the  line  and  region  files. 

Name  -  'S'  or  *L'  with  automatically  assigned 

number; 

Window  with  automatically  defined  dimensions; 

Indices  (addresses)  of  start-  and  endpoint  in  the 
coordinate  list; 

Coordinates  of  start-  and  endpoint. 

For  quality  control  and  to  support  the  2nd  step  the  data 
can  be  output  in  the  form  of  lists  or  graphical  plots  (pro¬ 
gram  MAPOUT,  see  section  3.3.2). 


(b)  Step  2:  Graph  generation 


Program  GRAPHREAD  generates  the  graph,  fills  up  the  region 
file  by  adding  to  it  all  areas,  and  completes  the  line-  and 
region  files  with  all  pointers. 

Input  to  program  GRAPHREAD  at  this  time  is  manual.  It 
consists  of  names  and  types  of  regions  and  identification  of 
area  boundaries  as  a  sequence  of  numbers  of  lines  and  pseudo¬ 
lines.  Limited  resources  prevented  us  from  going  into  auto¬ 
mation  of  graph  generation. 


(c)  Printout  for  control 


The  file  sizes  can  be  easily  verified  in  a  specific 
control  printout.  Figure  3.9  presents  a  section  of  such  a 
control  printout  of  the  data  bank  GRAZ.  Units  are  computer 
words  (36  bits)  as  occupied  on  disk.  The  logical  unit  of  the 
coordinate  file  is  a  pair  of  coordinate  values  occupying  2 
words .Therefore  there  are  1910  coordinate  pairs  in  the  file. 

The  logical  uniut  of  the  region-  and  line  files  in  28  words.  There 
are  thus  6216:28  =  222  regions  and  3304:28  =  118  lines. 
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Figure  3.9:  Control  printout  of  map  data  bank  Graz. 
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3.3.2  Data  output 


Data  output  is  with  program  MAPOUT  in  the  form  of  lists 
and  graphical  plots.  Output  routines  are  interactively  star¬ 
ted  (format  free  input  with  input  requests) .  Output  devices 
available  at  this  time  are  line  printers  and  drum  or  flat  bed 
plotters . 

Program  MAPOUT  was  used  to  generate  figures  3.2  through 
3.7  and  3.9.  It  can  also  produce  combinations  of  lists,  such 
as  of  coordinate  and  line  lists  as  shown  in  figure  3.10. 
Graphical  plots  in  figure  3.3(a)  were  enlargements  of  speci¬ 
fied  features  of  the  data  bank  of  which  a  larger  portion 
was  presented  as  figure  3.2,  all  plotted  with  the  same  pro¬ 
gram  MAPOUT. 


3.4  CONCLUSION 


Several  structures  exist  for  the  digital  map  data  in  the 
current  application.  These  can  be  grouped  into  three  main 
groups : 

-  the  map  data  bank  with  its  coordinate  files, 
element  lists  and  topological  structure; 

-  the  polygon  (vector)  structure  of  individual  map 
features ; 

-  the  (binary)  matrix  structure  for  compatibility 
with  digital  images. 

The  apparent  complexity  of  a  digitally  implemented  topological 
structure  has  great  benefits  when  the  data  must  be  analyzed 
and  lengthy  search  processes  can  be  drastically  shortened. 

The  polygon  (vector)  structure  of  individual  map  features 
is  an  intermediate  structure  that  is  feature-oriented.  It 
has  the  drawback  of  wasted  storage  if  large  data  quantities 
need  to  be  kept  in  this  format.  But  in  otxr  application  this 
drawback  is  irrelevant  because  the  structure  is  only  used  as 
an  intermediate  one  for  the  ultimate  presentation  in  a 
binary  image;  this  is  the  structure  needed  to  input  map  data 
into  a  digital  image  processing  system. 
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4.  MAP  DATA  IN  IMAGE  COMPATIBLE  FORM 


4.1  GENERAL 


The  main  current  application  of  the  map  data  bank  of 
chapter  3  is  to  serve  for  experiments  to  automatically  merge 
an  image  with  a  map.  Applications  programs  therefore  address 
the  problem  of  extracting  data  from  the  data  bank  in  a  format 
suitable  to  use  with  digital  images. 

Within  the  overall  flow  of  operation  this  task  occurs 
in  program  MAPTRANSFORM  which  in  turn  is  called  in  statement 
(9)  of  procedure  OVERLAY  (figure  2.5). 

The  application  programs  can  be  grouped  according  to 
the  data  structure  used,  figure  4.1  presents  an  overview. 


RAND  TRAPOL  POLBIN  ADDBIN 


NOTBIN 

SUBBIN 


Figure  4.1:  Existing  application  programs  for  the  map  data  bank. 


Programs  RAND  and  POLBIN  serve  to  convert  one  data  structure 
into  another  one.  The  other  programs  are  to  operate  on  a  given 
data  structure:  TRAPOL  performs  a  geometrical  transformation, 
set  operations  are  feasible  with  programs  MULBIN,  ADDBIN, 
NOTBIN,  SUBBIN?  program  MASS  computes  the  centroid  of  a  feature 
from  either  its  coordinate  list  or  its  binary  matrix.  One  may 
argue  that  MASS  creates  a  new  data  structure  (the  center  of 
gravity)  from  other  structures . 

Program  RAND  extracts  from  the  data  bank  the  contour  of  a 
feature  or  object  in  the  form  of  a  coordinate  list.  These  da¬ 
ta  are  then  available  in  either  a  vector  format  or  in  raster 
presentation  for  further  processing  in  the  context  of  digital 
image  analysis . 


4.2  EXTRACTION  OF  A  FEATURE  FROM  THE  DATA  BANK/ PROGRAMM  RAND 

As  was  mentioned  before  it  is  through  routine  MAPTRANSFORM 
that  a  map  feature  is  transformed  into  an  image  compatible  form. 
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The  PDL  formulation  is  presented  in  figure  4.2.  It  calls 
the  programs  listed  in  figure  4.1. 


(1) procedure  MAPTRANSFORM  (map,  object,  binarymatrix); 

(2)  begin 

(3)  RAND  (map,  object,  polygon  1); 

(4)  TRAPOL  (polygon  1,  polygon  2); 

(5)  POLBIN  (polygon  2,  binarymatrix) ; 

(6)  end 


Figure  4.2:  PDL  formulation  of  program  MAPTRANSFORM 


Program  RAND  is  important  for  ARSIM  because  it  addresses 
individual  features  of  the  map  data  bank.  Given  the  name 
Cregnam)  of  an  object  RAND  generates  the  boundary  polygon 
as  a  stream  of  x,y  coordinates. 

The  PDL-description  of  RAND  is  presented  in  figure  4.3. 
RAND  first  uses  program  KASUCH  (map,  regnara)  to  search  in 
the  region  file  for  the  number  NREG  of  the  description  of  ‘ 
the  object  and  from  there  takes  the  REG-item  associated  with 
object  REGNAM.  Should  the  description  of  the  region  show 
that  it  is  a  skeleton,  then  the  search  is  simple;  a  file  is 
filled  with  a  sequence  of  skeleton  points  from  the  skeleton 
coordinate  list  (Fig. 3.1)  in  positive  and  negative  direction 
(printl(+l),  printl(-l)).  This  ensures  that  in  a  conversion 
to  a  binary  matrix  one  has  2  border  points  along  each  direc  - 
tion  of  a  raster  line.  This  also  enables  one  later  to  widen 
the  skeleton. 

Is  the  object  an  area  then  RAND  goes  along  all  LIST-  and 
LIN-items  of  the  region  using  pointer  IREG  and  routine 
CONTINUOUSLINE.  The  main  subroutine  in  RAND  is  CONTINUOUSLINE . 
It  is  presented  in  figure  4.4. 

In  order  to  describe  RAND  and  CONTINUOUSLINE,  let  us 
now  consider  two  different  cases,  first  a  simple  one,  with  an 
area  that  is  only  defined  by  lines,  then  the  complex  one, 
with  an  area  within  another  area,  defined  by  lines  and 
skeletons . 


4.2.1  Area  defined  only  by  lines 

The  search  in  the  region  file  results  in  a  pointer  IPL 
(=  element  (c)  in  chapter  3.2.6)  to  a  REG-itm  (IREG).  This 
in  turn  produces  the  first  LIST-item  ('liste')  and  enters  into 
subroutine  CONTINUOUSLINE.  This  defines  through  routine  REGLIS, 
whether  or  not  the  associated  LIST/LlN-item  concerns  a  line 
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(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 
(11) 
(12) 

(13) 

(14) 

(15) 

(16) 
Cl  7) 
(18) 

(19) 

(20) 
(21) 
(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


procedure  RAND  (map,  regnam,  polygon) ; 
begin 

LPS0:=0/ 

LISTA: =0/ 

NREG :  =KASUCH  (map,  regnam); 
if  NREG  <0  then  "not  found" 
else 

begin 

readrg  (SI,  NREG); 

IREG: =IPL  (SI)/ 

if  RA  (IREG)  NREG  then  "DB-error" 

else 

begin 

assign  polygon  file; 
initiate  polygon; 
if  regnam  is  a  skeleton  then 
begin 

printl  (+1); 
printl  (-1) ; 
close  polygon; 
close  polygon  file; 
end  else 

begin 

liste:=  PLIST (IREG) ; 
while  liste  4  0  do 
begin 

CONTINUOUSLINE  (LPS0,  LISTA,  liste); 
close  polygon; 
liste ;=LISTA; 

LPS0 : =LISTA : =0 

end 

close  polygon  file; 

end 

end 

end 

end. 


Figure  4.3;  PDL  description  of  routine  RAND. 


(1) 

procedure  CONTINUOUSLINE  (LPS0.LISTA, liste) ; 

(2) 

begin 

(3) 

case 

REGLIS  (liste)  of 

(4) 

'L* : 

begin 

(5) 

readl(Sl) ; 

(6) 

case  REGLIS(liste)  of 

(7) 

'L' :  begin 

(8) 

readl(S2) ; 

(9) 

switch:*  1; 

do) 

if  start  l*start  2  or 

(ii) 

begin 

(12) 

printl(-l) ;top:*en 

(13) 

end 

(14) 

begin 

(15) 

top:*start  1;  prin 

(16) 

end 

(17) 

end 

(18) 

1 P'  switch:*  2; 

(19) 

eof:begin  printl(+l) '.switch:*  3 

(2o) 

end 

(21) 

'P*  . 

begin 

(22) 

readrg  (SI); 

(23) 

LPS0:  ■  line; 

(24) 

case  REGLIS  (liste)  of 

(25) 

'L':  begin 

(26) 

readl  (S2); 

(27) 

switch:  *  1; 

(28) 

top  1:  ■  top:  *  NAEHER 

(29) 

end 

(3o) 

' P' :  switch:  *  4; 

(31) 

eof:  begin  printl(+l) ;  switch:* 

(32) 

end 

(33) 

eof : 

begin 

(34) 

if  LISTA  ■  0  then  "no  boundary" 

(35) 

switch:  ■  3; 

(36) 

end 

(37) 

case 

switch  of 

(38) 

1:  begin 

(39) 

if  top  ■  end  2  then  printl(-l) 

(4o) 

else  printl(+l); 

(41) 

L4; 

(42) 

end 

(43) 

2:  begin 

(44) 

readrg  (S2); 

(45) 

top:*  NAEHER  (start  l.end  1,S2); 

(46) 

if  top*start  1  then  printl(-l) 

(47) 

else  printl(+l); 

(48) 

L5; 

(49) 

end 

else 


Figure  4.4.  Routine  CONTINUOUSLINE(cont 'd  on  nextpage). 
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(50) 

(51) 

(52) 

(53) 

(54) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

(60) 

(61) 
(62) 

(63) 

(64) 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 

(71) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 

(77) 

(78) 

(79) 

(80) 

(81) 
(82) 

(83) 

(84) 

(85) 

(86) 

(87) 

(88) 

(89) 

(90) 

(91) 

(92) 

(93) 

(94) 

(95) 

(96) 

(97) 

(98) 

(99) 

(100) 

(101) 

(102) 


3:; 

4:  begin 

if  IREGA-IREGB  Chen 
begin 

case  SPITZ  of 
1:  begin 

top  l:*top:“start  1; 

S2:-S1; 

end 

2 :  begin 

top  l:*top:“end  1; 

S2:-S1; 

end 

3:  begin 

LPS0:  -  0; 

case  REGLIS  (liste)  of 
'L' :  begin 

readl  (S2); 

if  start  l“start2  or  start  l“end  2  then 
begin 

print 1 (+1 ); print 1(-1) ; topi: “start  1 
end  else 

begin 

top  1 :*endl ;  printl(-l);printl(+l) 

end ; 

if  top“end  2  then  print 1(-1) 
else  printl(+l); 

1*4; 

switch :*1;0  dont  enter  L5  in  line  loo  % 

end 

'P' :  begin 

readrg  (S2); 

if  start  l“start  2  or  start  l*end  2  then 
begin 

print 1(+1); printl (-1 ); topi : “start  1 
end  else 

begin 

top  1 :*endl ; print 1  (-1) ;printl(+l) 

end; 

IFEGA: “IREGB 

end 

eof:  begin  "Ireg“skeleton"; switch:*  3  end; 

end 

end  else 

begin 

t  IREGA-IREGB  0 

readrg  ( S2 ) ; 
top:“start  1; 

if  top-start  2  &  top-end  2  then  top:*end  1; 
top  1:«  top; 

IREGA:  “IREGB 
end 

if  switch  ■  4  then  L5; 

end 
end  • 


Figure  4.4.  Routine  CONTINUOUSLINE(contil»ued) . 


(1)  procedure  L4; 

(2)  begin 

while  REGLXS  (liste)  =  'L'  do 

(4)  begin 

(5)  readl  (Si) ; 

(6)  if  top  *  endl  then  printl  (-1) 

(7)  else  printl  (+1); 

( 8 )  end 

(9)  if  REGLIS  (liste)  =  'P'  then 

(lo)  begin 

(110  .  readrg(s2); 

(12)  L5 

(13)  end 

(14)  else 

(15)  if  REGLIS  (liste)  =  eof  then  L6; 

(16)  else  error; 

(17)  end. 


Figure  4.5:  PDL-description  of  procedure  L4  within 
CONTINUOUSLINE .  L4  follows  a  contour 
along  lines  (segments) . 
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( '  L ' ) ,  a  pseudo-line  ( ' P ' )  or  the  end  of  the  feature  C ' eof ' ) . 
With  the  present  assumptions  it  must  be  a  line  'L'.  We  have 
to  find  now  the  direction  in  which  the  sequence  of  coordinates 
must  be  ordered. 

Therefore  one  reads  from  the  line  list  the  current 
line  (readl  (SI))  and  the  subsequent  one  (readl  (S2));  SI  and 
S2  denote  two  different  storage  buffers.  Then  it  is  a  trivial 
task  to  find  the  direction  from  statement  (10)  of  program 
CONTINUOUSLINE .  One  now  has  to  copy  the  points  of  the  line  from 
the  coordinate  file  to  the  polygon  file  in  the  order  resulting 
from  statement  (10),  either  as  printl  (+1)  or  printl  (-1); 

"top"  is  set  as  the  last  point  read  and  is  kept  for  further 
use.  successive  lines  are  then  added  on  using  procedure  L4 
(f'igure  4.5).  This  is  repeated  until  pointer  PNEXT  of  the 
LIST-item  "liste"  equals  zero.  Then  procedure  L6  is  executed 
(f'igure  4.8)  to  terminate  the  polygon. 


4.2.2  An  area  within  another  area  defined  by  lines  and 
skeletons . 


This  is  the  more  complex  case.  It  becomes  relevant  if  a 
LIN-item  represents  a  pseudo-line.  One  of  the  two  pointers 
PR1,  PR2  must  then  point  to  the  parental  REG-item,  the 
other  pointer  can  also  point  there.  In  that  case  one  is 
dealing  with  a  pseudo-line  connecting  for  example  one  island 
with  the  boundary  of  a  lake  (area  within  area) .  Or  the 
other  pointer  must  point  to  a  skeleton.  Pseudo-lines  connect¬ 
ing  islands  with  a  region  boundary  are  detected  in  proce¬ 
dure  REGLIS.  Since  they  do  not  affect  the  polygon  of  the 
current  boundary  level  the  part  of  the  line  list  between 
the  first  and  the  second  appearence  of  such  pseudo-lines 
are  skipped  in  REGLIS  and  the  type  of  the  following  LIN- 
item  is  returned.  To  resume  the  skipped  list  part  the 
pointer  to  the  first  LIST-item  after  the  pseudo-line  is 
retained  in  LISTA.  After  completion  of  the  current  polygon 
the  polygon  circumscribing  the  island  is  added  to  the  poly¬ 
gon  file  as  shown  in  lines  (25'*  to  (31)  of  procedure  RAND. 

If  a  feature's  contour  meets  with  a  skeleton  the  al¬ 
gorithm  enters  procedure  L5  (Figure  4.6).  The  entry  state 
is  defined  by  the  "top"  of  the  contour  generated  so  far 
and  a  skeleton  stored  in  buffer  S2.  The  direction  in  which 
the  contour  follows  the  skeleton  is  determined  by  the  next 
LIST-item.  There  are  3  cases  to  distinguish. 


Case  SKI  (figure  4.7): 

From  the  skeleton  border  the  contour  continues  with  a 
line  (lines  49-55  of  L5) .  Rename  "top  0".  Procedure 
NAEHER  determines  the  new  top  which  is  either  the  start¬ 
ing  (start  1)  or  the  end  point  (end  1)  of  the  line  so 
that  it  is  nearer  to  the  skeleton.  Then  the  polygon  along 
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(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(10) 

(11) 
(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 
(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


procedure  L5? 
begin 

LIND:=  true; 

while  REGLIS (liste) =  *P'  &  LIND  do 
begin 

if  I REGA= I REGB  then 
begin 

case  ISPITZ  of 
1 : begin 

RASKEL (top, start  2,S2); 
top:=start  2 

end 

2: begin 

RASKEL (top, end2 ,  S2 )  ; 
top:=end  2 

end 

3 : begin 

if  top  =start  2  then 
begin 

.  RASKEL (top, end2 , S2) ; 
RASKEL (end2, top, S2) ; 
end  else 

begin 

RASKEL (top, start2,S2) ? 
RASKEL ( start 2 , top ,  S2 )  ; 

end; 

LIND;=false 

end 

end  else 


Figure  4.6.*PDL-description  of  procedure  L5  within 

CONTINUOUSLINE.  L5  follows  a  contour  along 
pseudolines , which  separate  the  region  from 
adjacent  skeletons  (continued  on  next  page) . 
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(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 


begin  %  IREGA  *  IREGB  $ 
readrg (SI) ; 

if  start  1=  start  2  or  start  2=end  1  then 
begin 

RASKEL (top, start2,S2) ; 
top:=start  2 

end  else 

begin 

RASKEL (top , end2 , S2 )  ; 
top:=end2 

end 

IREGA:  =  IREGB; 

S2 : =Sl 

end 

end  while; 

if  not  LIND  then  L4 

else 

if  REGLIS  (liste) =  'L'  then 
begin 

readl  (Si) ; 
top  0:  =  top; 

top :=NAEHER( start  1,  endl,  S2)  ; 

RASKEL ( top  0,  top,  S2); 
if  top=start  1  then  printl(+l) 
else  printl(-l); 

L4 

end  else 

if  REGLIS (liste)  =  eof  then  L7 

else  error 


end 


Figure  4.6:  PDL-description  of  procedure  L5  within  CONTINOUS- 
LINE.  L5  follows  a  contour  along  pseudolines, 
which  separate  the  region  from  adjacent  skeletons 
(continued  from  previous  page) . 


skeleton 


contour 


Figure  4.7:  Case  SKI  of  a  skeleton  meeting  a  contour. 


the  skeleton  between  top  0  and  top  is  generated  by  proce¬ 
dure  RASKEL.  From  top  the  algorithm  runs  as  described  in 
chapter  4.2.1  (procedure  L4) . 


(1)  procedure  L6; 

(2)  begin 

(3)  if  LPS0  7*  0  then 

(4)  begin 

(5)  IREGA:  *  PR1  (LPS0); 

(6)  if  IREG  =  IREGA  then  IREGA  =  PR2  (LPS0) ; 

(7)  readrg  (SI); 

(8)  RASKEL(top, topi ,  SI)  ; 

(9)  end; 

(10)  end* 


Figure  4.8: 

i 


PDL-description  of  procedure  L6  within  CONTINUOUS 
LINE.  L6  closes  a  contour  when  coming  from  L4. 


contour 


/ 


Figure  4.9: Case  SK2  of  a  skeleton  meeting  a  contour 


Case  SK2  (figure  4.9): 

It  follows  a  pseudo  line  of  the  same  skeleton  CIREGA=» 
IREGB;  lines  6-29  of  L5) .  This  happens  for  example  when 
a  dead  end  road  enters  the  region.  The  contour  must  fol¬ 
low  the  skeleton  until  its  limit  and  return  on  the  other 
side  of  the  skeleton.  The  two  limits  of  a  skeleton 
are  marked  in  the  LIST  by  special  PL-pointers:  262140 
marks  a  start  and  262141  an  endpoint.  Function  ISPITZ 
determines  which  limit  enters  the  region.  It  returns  1 
if  it  is  a  start  point,  2  if  it  is  an  end  point  and  3  in 
a  special  case:  The  contour  joins  the  skeleton  in  a  limit 
point  and  the  other  limit  lies  inside  the  region.  This 
skeleton  is  described  by  exactly  2  pseudolines  and  the 
two  limit  marks .  Therefore  the  contour  must  follow  the 
skeleton  to  the  opposite  limit  and  return  back  to  the 
starting  point,  which  remains  the  top  of  the  contour 
polygon.  Hence  the  algorithm  returns  to  procedure  L4. 

In  the  first  two  cases  ISPITZ  =  1  and  ISPITZ  =  2  the  way 
of  the  contour  can  be  traced  until  the  limit  point,  which 
becomes  the  new  top.  The  algorithm  is  then  in  the  same 
state  than  before  and  iterates  in  procedure  L5 . 
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(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(lo) 
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(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 
(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


procedure  L7; 
begin 

if  LPS0?0 then  RASKEL (top. top  1,S2) 
else 

begin 

IREGA=PR1  (LPS0)  ; 

if  IREG=IREGA  then  IREGA:=PR2 (LPS0) ; 
if  IREGA=IREGB  then 
case  ISPITZ  of 
1 : beg in 

RASKEL (top, start2,S2)  ; 

RASKEL (start 2, top  1,S2) 

end 

2 : begin 

RASKEL (top, end2,S2)  ; 

RASKEL (end2, top  1,S2) 

end 

3 : if  top=start2  then 
begin 

RASKEL ( top , end2 , S2 ) ; 

RASKEL (end2, top  1,S2) 
end  else 

begin 

RASKEL ( top , start  2 ,  S  2 )  ; 

RASKEL (start2, top  1,S2) 

end 

else 

begin  0  IREGA^IREGB# ; 
readrg  (SI); 
top0:»start  1; 

if  end  l=start  1  or  end  1=  end  2  then 
top0:=end  1; 

RASKEL (top, top0,S2) ; 

RASKEL (top0, top  1,S1) 

end 

end 

end- 


Figure  4.10;  PDL-discription  of  procedure  L7  within  CONTINUOUS 
LINE.  l7  closes  a  contour  when  coming  from  L5. 
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Figure  4.11s  Case  SK3  of  a  skeleton  meeting  a  contour. 


Case  SK3  (figure  4.11): 

It  follows  a  pseudo-line  of  another  skeleton  ( IREGA=IREGB ? 
lines  30  -  43  of  L5) .  The  contour  then  follows  the  first 
skeleton  (by  RASKEL)  until  the  joint  point  with  the  other 
one,  which  is  the  new  top  to  continue  and  interate  in  L5 . 

There  exists  a  case  where  the  end  of  the  LIST  follows  the 
skeleton;  this  is  not  treated  here  separately,  because  a 
contour  is  always  a  closed  line.  The  first  LIST-item  there¬ 
fore  follows  the  last  one  and  allows  to  distinguish  the 
three  cases.  Some  detail  that  is  different  here  from  the 
previous  one  is  shown  in  procedure  L7  (figure  4.10). 


4.2.3  Calling  program  RAND 


Figure  4.12  is  an  example  of  an  interactive  call  to 
program  RAND  with  the  response  as  an  operator  receives 
it.  Two  areas  are  called: 

EGG2NBERG  and  SCHLOSSEGGEN 


that  we  have  used  in  chapter  3. 


J 
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•«  scon 

BEGRFNZUNGSPOL YfiONE  BESTINhEN 

«*  S  BAND  iiRA 2 •  EGGE N8E  RG 

FEN  5TE 9  WON  EGGENBE  PG  = 

K00R0IN*T£MFIL£  EGGENBERG 

21 1. 746  20*6 .027 

FRSTELLT. 

5.**  Jo 

9.997 

elapsed: 

3.oro  TOTAl  : 

2.5  K4  CHt:  .  1J.5 

i/o: 

1.828  SEC 

**  iPANH  SCHlOSSEGl.E* 

FENSTER  VON  SCHl OSSt GGF N= 

kooroinatenfile  schlosseggfv 

,?13.6  9ff  2C  48.16“ 
ERSTELLT. 

1. 224 

2.12n 

elapsed: 

.773  total: 

1.859  CPU!  .02° 

i/o: 

1.202  SEC 

Figure  4.12:  Call  to  extract  regions  EGGENBERG  and  SCHLOSS- 
EGGEN  from  the  data  bank  GRAZ. 


The  call  includes  the  name  of  the  data  bank  (GRAZ)  and 
of  the  region  (e.g.  EGGENBERG) .Interactive  commands  are 
preceded  by  a  sign.  The  printout  also  generates  in**'. 

The  machine  response  has  no  preceding  symbol.  Times  given 
in  the  last  row  of  figure  4.12  are  machine  dependent.  On  our 
machine  (UNIVAC  1100/81)  the  CPU-time  was  0.029  seconds. 


4.3  OPERATIONS  WITH  THE  EXTRACTED  FEATURE 
4.3.1  TRAPOL 


Procedure  TRAPOL  (figure  4.13)  performs  a  simple  linear 
geometric  transformation  of  the  polygon  created  by  program 
RAND.  It  also  can  widen  individual  skeletons  to  a  specified 
width . 

The  transformation  T(x)  is: 

1(X)  sAl  +  i 

with  vectors  Ks),  £and  b  C  R2  and  a  square  transformation 
matrix  &.  Transformation  parameters  must  be  computed  in  a 
separate  program. 
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(1)  procedure  TRAPOL  (polygon  1,  polygon  2) ; 

(2)  begin 

(3)  read  parameters  of  transformation  T; 

(4)  assign  file  for  polygon  2; 

(5)  initiate  file; 

(6)  for  all  part  of  polygon  1  do 

(7)  begin 

(8)  for  all  point  £  of  part  do 

(9)  polygon  2;  =  polygon  2  +  T  (jg) 

(10)  close  polygon; 

(11)  end 

(12)  close  file; 

(13)  end. 


Remarks  to  procedure  TRAPOL: 


1)  Description  of  a  polygon  as  a  syntagmatic  grammar 
(VAUQUOIS,  1970): 

2 

G  -  (  1R  VJ  {  endofpart),  (polygon,  part,  point)  ,  polygon,  P  ) 
with 

P  =  (polygon  —-♦part  polygon  /  part 

part  — — ♦  point  part  endofpart  /endofpart 

point  — ►  (x,  y)  £lR2} 

2)  "polygon  +  point"  means  attaching  the  point. 

3)  "close  polygon"  writes  the  endofpart. 

4)  T  (jg)  =  +  b  with 


.  /aU  a12  \ 
\  a21  a22  / 

•  (y 


and 


Figure  4.13:  Procedure  TRAPOL 
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Widening  a  skeleton  is  achieved  by  displacing  each  part 
of  the  polygon  in  a  positive  direction  for  one  half  of  the 
required  width.  Polygon  elements  are  then  intersected  and 
the  intersection  points  are  then  used  as  the  new  polygon  points , 
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206.  2040. 
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1.  1. 

60  • 
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TX 
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A .2 1 A*X  ♦ 
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- 
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-2.341*Y 

4 

A8 16.800 

•* 
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8EGREN7UN6SP0LYG0 W  INS 

NEUE  K00R0INATE  NSYSTEH  t  kansforhieren  : 

STRAPOL  EGGENBERG* 

EGPOL 

transforhierter  file 

EGPOL 

HAT 

168  PUNKTE 

SEIN 

fenster  : 

25  3 

25 

27 

elapsed: 

•  55A  total:  1.083 

cpu: 

.015  I/O:  .453  SEC 

Figure  4.14:  Call  of  linear  transformation  of  region  EGGENBERG. 


»*  STRapol  Sr.  jOCt.ctCCAS*  S45FPL 
TRA»'bf  TRNIFPTER  FILE  S45P0L  HAT  93  HUwKTt 

SEIM  FENST FR  :  2  24  41  IQ 


ELAPSE o: 


.  a  to r : 


.952  CPU: 


.f’16  i/o: 


326  SFC 


**  scow  SKELETON  VERB RE  ITf  Rrw  <  STPASSE  2  X  0.3  ) 


**  *r*APOL  s  *»  ra  p  I  -  L  tfA  ^PjL  -4  #  .3 
TRAMSFORNIERTEF  FILF  S45P01-4  HAT  1C2  HUNKTE 

SEIM  FENSTFP  :  ?  24  41  10 


FLA^SF.0: 


.TPA  TOTAI  I 


.(iso  cpo: 


."2f  i/o: 


.227  iiC 


Figure  4.15:  Call  of  program  TRAPOL  to  transform  and  widen 
the  skeleton. 


L . A 
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Figure  4.14  is  an  example  on  an  interactive  call  to 
the  linear  transformation.  Here,  transformation  parame¬ 
ters  are  computed  using  3  control  points  in  an  interac¬ 
tive  call  to  program  MASS  that  is  then  followed  by  TRAPOL. 
MASS  uses  3  points  (6  coordinates)  in  the  map  system, 
followed  by  3  points  (6  coordinates)  in  the  transformed 
system.  The  call  to  TRAPOL  requires  two  coordinate  lists: 
the  one  to  be  transformed  and  the  one  after  transformat¬ 
ion  (EGPOL) . 


Figure  4.15  is  an  example  for  transforming  skeleton 
S00000000045  into  S45POL,  and  then  widening  it  to  0.3 
units  and  storing  it  in  file  S45POL-4. 


4.3.2  POLBIN 


A  discussion  of  procedures  is  decribed  by  NEWMAN  and 
SPROULL  (1979) .  The  procedure  used  in  our  case  is  illustrat¬ 
ed  in  Figure  4.16.  One  intersects  each  part  of  the  polygon 
with  equidistant  (vertical)  straight  lines.  Intersection 
points  we  ordered  along  the  partial  polygon  (horizontally) , 
and  are  then  also  ordered  within  each  straight  (vertical) 
line.  One  then  has  a  raster  image  of  the  contour. 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 
(ID 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
(21) 
(22) 


procedure  POLBIN  (polygon,  binmatrix) 
begin 

binmatrix:  =  (0)  ? 

for  all  part  of  polygon  do 

begin 

binloc:  =  (0) ; 

P0:  =  first  point  of  part: 

for  all  point  PIC  part  -  |P0}  do 

begin 

if  x(Pl)  -  x(P0)  *  0  then 

for  i;=y(p0)  until  y(Pl)  do 
binloc  (x(P0),i):=  1; 

else 

INSORT  (store,  P0,  PI); 

P0:  »  PI? 

end 

for  all  columns  IX  of  binloc  do 
while  3  ya,  Ye  £  store  do 

for  iy:  =  [ya+.5]  until  [ye+.5]  do 

binloc  (IX,  IY) :  =  1? 
binmatrix:  *  |  binmatrix  -  binloc  J  ? 

end 
end  . 


Figure  4.16:  Procedure  POLBIN. 
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It  is  possible  to  number  all  edge  raster  points  on  a 
straight  line  and  set  all  raster  points  between  odd-numbered 
and  even-numbered  given  edge  raster  points  to  1,  all  others 
to  0.  This  then  gives  a  binary  matrix  of  the  area,  a  so-called 
"fill-in" . 


**  *  COM  BE  GRENZUNG  SPOL  YGONF  IN  B IN  A£RE  MAjklZEN  UMFOKMEN: 


**  SP0LBIN 
A  US  EGP0L 

EGPOL  fEGBINtlU 
EGBIN 

MIT  FENSTER 

1 

1 

56 

75  GENERIEkT 

ELAPSEO: 

.992  total: 

1.075  cpu: 

.061 

i/u: 

.391 

SEC 

Figure  4.17:  Call  to  program  POLBIN  to  convert  vector  to 
raster  data . 


Figure  4.17  is  a  sample  call  to  program  POLBIN.  The 
binary  matrix  is  to  be  stored  in  file  EGBIN.  IW  is  a 
window  of  56  by  35  pixels.  If  IW  is  not  given,  then  POLBIN 
itself  computes  the  minimal  window  to  fit  the  binary  ima¬ 
ge.  The  binary  result  is  shown  in  Figure  4.18. 


SPM.TEN  0-01  ais  005* 


>••••••••••••••  ••  •  •  ••  ••  •< 

. . . XX 


0001  . 

oooz  . . 

0003  . . 

0004  ........ 

0005  ........ 

0006 

0007  . . .  XXX XXX X. 


. . 

... xxxxx. xxx. 
.xxxxxxxxxxxxx.  ..  ..  ..  .. 


0008 

0009 

0010 

9C11 

0012 

0013 

C01* 

0015 


C01 6  . . . .  ..X  XXXXXXXXXXXXX  XX  XX 


X  XXX  XX  ..  ..  —  . 

X XXX  XXX.  ..  ..  .. 
.XXX XX  XX..  ..  . 
..XXXXXXX.  M 
..  .X  XX  XX  ..  —  .■ 
....  XX  XX  ..  *.  .i 
....  .X XX XX X.  .i 
..  .XXX XXXXXXX. 


0017 

0018 

0019 

9020 

0021 

0022 

0023 

902* 


.  ..xxxxxx. 

iXXXXXX..  . 

..xxxxxx 

.  .XXXXXXX. 

.  ..xxxxx. ., 


h  .XXXXXXX*.  ..  m  .. 'm 


—  xx  xx  xxxxx. 

—  xxxxxxxxx. 

—  xxxxxxxx.. 

—  xxxx xxxxxx. 
..xxxxxxxxx.  xxxx  xxxxxx. 
.xxxxxxxxxxx xxxxx. . 
xxxx  xxxxxx  xx  xxxxx.  < 

0025  . .....XXXXXXXXXXXXXXXXX.  ..  , 

0026  . . XXXXXXXXXXXXX!  ..  ..  -  .. 

0027  . . XXXX  XXXX  XX.  -  , 

0028  ...................... ....XXXXXX....  —  ....  ..  .. 

0029  . .  . . . 

0030  .....................................  »*..«•' 


i 

i 

I 

I 

I 


Figure  4.18:  Example  of  binary  matrix  EGGENBERG 
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4.3.3  Geometric  set-operations 


The  raster  format  is  well  suited  for  the  most  important 
geometric  set-operations.  Two  or  more  binary  images  can  be 
connected  on  a  bit  level  using  logical  functions  that  are 
very  fast:  the  currently  used  computer  has  36  bits  per  word  and 
36  logical  operations  can  be  done  in  parallel.  Program 
MULBIN  forms  the  intersection  of  two  areas  (logical  operator 
AND),  ADDBIN  forms  the  union  (logical  operator  OR),  NOTHIN 
forms  the  negative  of  an  area  (logical  operator  NOT)  and 
SUBBIN  produces  the  absolute  difference  between  two  areas 
(EXCLUSIVE  OR) . 

Figure  4.19  shows  the  printout  of  the  program  call  and 
the  result  of  a  union  of  region  EGGENBERG  in  binary  format 
(stored  as  SCHBIN)  and  a  road  (S45BIN) .  Programs  LIBIN  and 
STOBIN  are  input/output  routines. 


4.3.4  Map  information  as  an  image 


The  binary  matrices  now  enable  one  to  generate  chi or o- 
plethmaps  of  the  regions  in  the  data  bank.  Be  M,  ,  i  =  1 , 

...,  n  mutually  exclusive  binary  matrices  and  g.,  i  =  1, 

...»  n  [n  G  0,255]  associated  "gray  values"  we  xthen  can 
build  a  gray- tone  digital  image  B: 

n 

B  -  £(g±  .  M±) 

i=l 


Image  H  can  now  be  processed  in  an  image  processing  system 
such  as  our  DIBAG  (LEBERL,  KROPATSCH.,  1979)  . 

An  example  of  such  a  chloropleth  presentation  is  shown 
in  figure  4.20  with  4  regions: 

areas  EGGENBERG,  SCHLOSSEGGEN,  HAHNHOF,  and  skeletons 
SOOOOOOOOO88,  S00000000045 . 


4.3.5  Program  MASS 


MASS  computes  the  center  of  gravity  S  of  a  region  and 
its  area  for  both  vector  or  raster  format.  The  formulae  are: 

(a)  vector  format  with  polygon  (x^  yi),  i  =  1,  ...  n 
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a) 

**  *COH  VEPEINIGUNG  VON  SCMLOSS  nIT  S4  5  : 


**  SLIBI*  SCHB IN 

BINAFRE  MATRIX  VON  FILE  SCH8IN  HIT  FENSTE*  1  l  56 


**  SAOOBIN  S 45 BIN 


**  SSTOBIN  VEREIMCUNG 

BINAERE  MATRIX  AUF  FILF  VE-VF  IN IGUmG  GESPEICHERT 


b) 

SPALTEw  Cf-Cl  BIS  00  56 


0001 
C002 
C003 
COC* 
0005 
0JC6 
0007 
000  H 
OoOP 
'••010 
.'Oil 
CO  12 
r'0t7 
uO  1 4 
<UU  5 
on  6 


CO  17  . . . •  •• . xxxxx.. 

'"•018  .  .... . .XXXXXX... 

<019  XXXXX... 

0020  XXXX... 

'-<•21  ....XXXXX... 

r022  . . . . •••••• . X... 

C023  . . . . . . 

ru24  ••••••••■•••••••••••••••••••••••••••••••• 

rnzs  ..•••••••••••••••••••.•••••••••••••••xxxx 

^026  . . . XX.. . XXXX.. 

f 027  .....................  XX xx XXXX. .XXXX...... 

•?t'2S  . . . . ...XXXX  XX  XX. ...... 

1*029  ..X.... ............ .XX. . . . . 

ng3C  ..XXX XXX. ....... ..XX.  ••••••  ..••••••..•••• 

...... xxxxxx. .. ..xx  ..  . . . 

•  mi  3?  ...•••••••x  xxxxxx  •••••••••.•••••••.  •••••• 

f.m  ....... ....xxxxx. ......  .................. 

r  M34  . . . 

•.-35  . . . . 


Figure  4.19.:  Call  to  program  ADDBIN  to  combine  area 
SCHBIN  with  skeleton  S45BIN  (a  road) , 
and  the  result  of  extracting  the  program. 
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Figure  4.20:  Chloropleth  map  of  Eggenberg. 


The  following  print  symbols  were  used: 


"0"  ...  Type  8,  railway  (BAHNHOF) 

"W"  ...  Type  7,  road  (skeletons) 

"O"  ...  Type  4,  open  area  (EGGENBERG) 

...  Type  3,  built  up  area  (SCHLOSSEGGEN) 
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area  (  (x^  yi)  ,  i  «  1#  . . .  n) 


2  E  (yi  •  Xi+1  -  yi+i  •  xi> 


i»i 


wlth  xn+l  5  X1  '  yn+l  5  yl 


(  (x^  y^)  *  i  —  I#  .  •  •  n) 


I  • 

3 


(y<  *x 


i--.:i-ti  "  yi-n,xi) 


n 

I 

i.l 


£  <i-i-xi*i  -  n-n^x' 


Sy  (  (xA>  Yj) •  1  =  1,  ...  n)  = 


n 

1  •  Ski  (x1'yH-l  -  xl+l-yl)  •  ) 

3  £  (xi-yi+l  -  xi+l-yi) 

i=l 

(b)  binary  matrix  B  (i,j) 

area  (B)  =  card  |B(i,j)  =  l| ; 

This  is  the  number  of  matrix  elements  of  B  equal  to  1. 


Sx  u> 


n  m 

(  £  H  (i  .  B  (i, j )  )  )  /  area  (B) 
i=l  J=1 


sy  caJ 


n  m 

(2  £  (j  •  B  (i, j )  )  )  /  area  (B) 

i-1  J-l 


Figure  4.  21  is  a  sample  call  to  program  MASS  to  illustrate 
its  result. 
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*•  *C0 N  SCHUtRPUNK  TE  ♦  FLAECHEN  VON  POLYGON  UNO  BI N  aft ATNI X  SCHLOSS  : 


•*  *L IB  I N  SCH8IN 

BINAERE  PA  TP  IX  VON  FILE  SCHB IN  NIT  FENSTEk  1  1  56  35  GEL 


**  *  NASS  SCHLOSSEGGEN*  BIN 

SCHLOSSEGGEN:  FL»ECHE=  1.577  schuerpunkt= 

BINAFRNATRIX:  FLA£CHE=  26.0C0  SC HUERP UNA T= 

TX  =  4.214*X  «■  • 330  *Y  ♦  -867.402 

TY  =  .000*X  ♦  -2. 341 *Y  ♦  4816.871 


214.421  2049.326 

36.231  19  .03  6 


**  *C0*  NEUER  VEPSCHIE8UNGSVEK  TOR  OURCH  GLEI CHSETZE  N  OEN  SCHwERPUNKTE 


Figure  .4.21:  Call  to  program  MASS  to  compute  area  and  center 
of  gravity,  and  corrects  the  transformation 
function. 


Binary  image  SCHBIN  is  first  read  in.  Then  the  call  to  MASS 
results  in  the  area  and  centers  of  gravity  of  the  vector 
presentation  (file  SCHLOSSEGGEN)  and  of  the  binary  presents  - 
tier  (file  SCHBIN) .  The  results  are  in  the  coordinate  system 
of  the  data  bank  (SCHLOSSEGGEN)  and  of  the  binary  matrix 
(SCHBIN) .  We  already  had  a  transformation  between  the  two 
shown  as  figure  4.13  with  3  manually  selected  control  points. 
The  center  of'  gravity  as  obtained  in  figure  4.18  is  now" 
automatically  transformed  and  if  a  discrepancy  with  the  ac¬ 
tual  center  of  gravity  in  the  other  data  struettire  exists  then 
a  shift  is  automatically  applied  to  the  transformation  (.com¬ 
pare  transformation  parameters  of  figures  4.13  and  figure  4.18). 


4.4  CONCLUSION 


It  has  been  shown  how  elements  of  a  map  data  bank  can 
be  extracted  and  converted  to  a  format  compatible  with  images: 
this  is  through  the  generation  of  the  contour  of  an  object 
from  the  data  bank  and  subsequent  conversion  to  a  binary 
image.  This  image  of  one  object  can  be  combined  with  those 
of  other  objects  to  a  digital  chloropeth  presentation  of  the 
map  contents.  This  can  then  be  processed  like  an  image: 
both  the  binary  and  the  chloropleth  presentation  can  be  un¬ 
derstood  as  synthetic  digital  images. 

The  binary  presentation  can  also  serve  as  a  mask  to  de¬ 
fine  a  specific  section  of  an  image.  This  application  is 
equally  important  then  that  of  serving  as  synthetic  image. 
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The  most  complex  searching  must  be  performed  in  program 
RAND.  However,  due  to  the  organisation  of  the  map  data  bank 
with  a  separate  planar  graph,  with  region-  and  line  files 
and  a  separate  coordinate  file,  searching  can  be  done  effi¬ 
ciently  and  fast. 
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5,  PRE-PROCESSING  OPERATIONS  WITH  IMAGES 


Any  form  of  digital  (or  analog)  image  analysis  bene¬ 
fits  from  an  optimization  of  the  image.  This  precedes 
the  analysis  and  is  denoted  by  "pre-processing" .  One  may 
classify  pre-processing  activities  as  follows: 

Single  images : 

Restoration 

Enhancement  and  Coding 

Geometric  Rectification 

Multiple  images: 

Data  Compression 

Enhancement  and  Coding  (e.g.  ratioing,  directional 

cosine,  clustering,  etc.) 

Of  the  countless  methods  of  image  (pre-) processing 
there  are  some  of  particular  significance  in  the  current 
context.  In  order  to  define  these  we  have  to  realize  that 
features  are  extracted  from  a  map  data  bank  mainly  in  the 
form  of  contoures  of  an  object?  these  are  then  presented 
in  a  binary  (black  and  white  only)  form.  It  follows  thus, 
that  also  the  multispectral  image  should  be  converted  to 
a  compatible  format.  This  can  be  achieved  by  a  clustering 
operation  of  image  pixels  (region  growing)  and  subsequent 
binary  coding,  to  enhance  specific  areas  or  the  contours 
and  edges. 

We  conclude  therefore  that  image  data  compression 
to  a  binary  form,  and  for  the  detection  of  edges  in  an 
image  ,  are  essential  pre-processing  steps  in  our  task. 

We  have  therefore  analyzed  a  number  of  methods  to 
identify  edges  in  images,  and  have  implemented  various 
data  compression  techniques  in  order  to  obtain  a  pre-pro¬ 
cessed  image  for  subsequent  line  detection  tasks.  This 
chapter  deals  first  with  a  mathematical  definition  of 
digital  images,  edges,  lines  and  neighbourhoods.  A  number 
of  pre-processing  functions  are  then  defined  in  mathemati¬ 
cal  terms  and  finally  evaluated  with  experimental  data. 


5.1  DEFINITIONS 


5.1.1  Digital  image 


We  describe  a  digital  image  by  a  function: 


q: 


B  -  G 


(5.1) 
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o£  a  finite,  k-dimensional  image  area  BCN^,  where  K  repre¬ 
sents  the  natural  numbers  1,2  ...  and  k  is  normally  2 . 

The  function  maps  B  into  a  finite  range  of  gray  values 
GCU.  For. multiple  images,  e.g.  multisprectral  ones, 
k  =  3.  We  can  define  multiple  images  also  conveniently  by 
a  vector  function 

g  -  g2'  •••  gn* 

g^  *  b  -► 

where  i  denotes  the  channels  of  a  multiple  image.  Each 
element  of  the  vector  function  can  be  represented  by  a 
matrix  of  gray  values.  These  matrices  consist  of  the  pixels 
defined  at  each  matrix  location  i  (row  ) ,  j  (  column  ) .  We 
may  in  the  following  denote  a  pixel  as  an  image  point  xCB 
or  by  its  location  ( i » j ) C  B . 


5.1.2  Neighbourhoods 

We  define  on  image  areaB  a  neighbourhood  I*(x)  to  each 
image  point  x. 

r(x)  C  Zk  is  the  set  of  neighbours  of  xGB  if 


x  £  T(x)  ,  (5.2) 

x  c  r(y)  y  c  r(x)  ,  (5.3) 


where  Z  is  the  set  of  integer  numbers. 

jrn(x)  C  B  j  n  £  is  a  system  of  neighbourhoods  around 
x  G  B  if: 


(a)  V  r(o)  (x)  :  =  {xf 
x  C  B 

(b)  r1  (x) :  *  r(x)r\B 

(c)  Tn+1(x):  = 


(5.4) 

(5.5) 


,^(.1  rly,nB 


(5.6) 


rn(x) 


It  is  to  be  noted  that 

rp“1(x).  r(x)  may  not  need  to  be  a  subset.of  B  wnile  rMx) 
is.r(x)  can  have  negative  numbers,  T  (x)  not. 


does  not  necessarily  contain 


There  exist  various  types  of  neighbourhoods  as  shown 
in  figure  5.2,  whereas  figure  5.1  illustrates  the  I\  neigh¬ 
bours  of  0th,  1st,  2nd  and  3rd  order  of  a  point  x. 


3 

3 

2 

3 

3 

2 

1,3 

2 

3 

_ 

3 

2 

1,3 

0*2 

1,3 

2 

3 

3 

2 

1,3 

2 

3 

3 

2 

3 

1 

3 

Figure  5.1:  ^-neighbourhoods  0tk,lst,2nc^  and  3r<^  order 
of  a  point  x. 


5.1.3  Metric 


We  need  for  further  use  a  metric  in  a  digital  image. 
All  pairs  of  points  (x,y)  ^B^  of  an  image  area  must  be 
related  through  a  neighbourhood  of  order  n.  A  distance 
between  x,y  is  then  defined.  We  can  prove  the  following 


Lemma  (1) : 

If  there  exists  for  each  point  x  C  B  a  system  of 
neighbourhoods  Ti  (x)  1  and 


y  a  b  =  LJ  ri(x)  (5.7) 

xCB  nCW  i=0 

O 


where  W  are  the  natural  numbers  and  zero,  then 
o 


d(x.y)  : 


min 


{kCBo 


x  e  rk(y)  | 


(5.8) 


is  a  metric  (distance)  on  B.  A  proof  can  be  found  in 
LEBERL  and  KROPATSCH  (1978b). 
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For  neighbourhoods  and  T.  of  Figure  5.2,  it  is 
easy  to  show  that 


d4(x,^):  =  £  /  xi  "  y±  / 

i=l 


(5.9) 


d8(x,^):  =  max  (/xi  -  y^/) 


(5.10) 


i=l , . .k 


B  1)  "4-Nachbar"  (Rosenfeld,  1970): 


Q  (  x  )  :  “(■]  Beachte:  card  (P4(x)  -  4 


B  2)  "8-Nachbar,r  KSnigszug  (Rosenfeld,  1970): 


m 


r8  C  x  )  •  Beachte:  card  (  Pg(  x  )  *  8 


B  3)  "Schrdge  Nachbaraengen"  konnten  bei  Bearbeitung  geometrischer 
verzerrter  Bilder  sinnvoll  sein: 


rT<x> .  -fij 


B  4)  Ala  Beispiel  einer  besonderen  Nachbarschaf t  seidie  "Springer"- 
beziehung  definiert: 


f  ( « )  <  -{el 


B  5)  Den  ZugmBglichkeiten  bain  "Dana"-Spiel  entspricht: 


rD  (  .  )  :  -  {■) 


Figure  5.2:  Five  examples  of  neighbourhoods  in  digital 
images. 
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5.1.4  Connectivity 

We  define  a  "path"  to  consist  of  image  points  p. , 
i  «1,...  n,  p^£  B  and  leading  from  p  ^  to  pn  if  1 

uL,  p*+icrl(p*>'“’lp*!  (5-U) 

We  define  a  subset  M  of  B,  M  C  B  to  be  "connected"  if 

V  x,y  CM3  pa^-h  (p^  from  x  to  yandVPi  C  M 

i 

and  a  subset  M  C  B  is  "simply  connected"  if 

-  M  is  connected;  _ 

-  M  =  B-M  is  connected,  whereby  in  M  a  different 
neighbourhood  relationship  may  be  defined  than  in  M. 

5.1.5  Edges,  lines,  regions 

We  now  have  the  tool  to  define  edges,  lines  and  regions 

(a)  An  "edge"  e  in  a  digital  image  is  a  pair  (i,o)  £  B2  of 
neighbouring  image  points  if 

i  €  T(o) .  (5.12) 

(b)  The  "weight"  or  "value"  of  an  edge  is  a  function 

2 

w:  B  R  CR  real  numbers).  (5.13) 

(c)  The  direction  of  an  edge  is  a  function 

OJ  B2  -*•  R  (5.14) 


We  see  that  each  pair  of  image  points  defines  an  edge. 
It  is  its  weight  that  provides  the  edge  with  a  physical 
meaning. 

One  needs  to  clearly  differentiate  between  edges  and 
lines.  We  therefore  define  a  "line"  L  as  a  finite  sequence 
of  edges  (e,),  i»l,  n  where  consecutive  edges  e,  ,  e,.,  are 
neighbours  of  one  another:  1  1+1 
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ei+l  €  A^ei^ 


(5.15) 


2 

where  symbol  A  denotes  a  neighbourhood  on  B  ,  analogous 
to  equations  (5.2)  and  5.3). 

We  can  define  open  and  closed  lines:  a  line  L=(e.,..e  ) 
is  'fclosed"  if 

el  £  A(en>- 
It  has  length  1  where 

l(L):»n  if  L  =  (e^  ,  i=l,  ...  n.  (5.16) 

A  "region"  in  a  digital  image  is  each  simply  connected 
subset  of  B. 

A  "contour"  L  of  a  region  R  is  a  line  L  =  (e^) ,  i=l,  ...n 
with  edge  e±  =  (p±,  q  )  £  B2  if 


V 


Pi  €  R  &  q±  ^  R 


(5.17) 


A  region  boundary  consists  of  image  points  b^.  i-=l.«,,n 
and  is  defined  if  for  all  edges  (p^,  q.)  =  e.  of  a  closed 
contour  L  =  (e^  ,  i=l,  n  of  region  R  one  finds: 


3  p  -  b.  4  J4i 
14j4n  1  2 


5.1.6  Comment  to  definitions 


It  is  significant  to  clearly  define  edges,  lines  and 
regions  to  avoid  misunderstandings  that  can  derive  from  com¬ 
mon  semantics.  Definitions  here  are  based  on  the  concept 
of  neighbourhood. 

A  possible  confusion  of  the  term  "edge"  may  commonly 
be  encountered  with  a  situation  as  shown  in  figure  5.3. 


ft 
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range  of  gray 
values  for 
R1  and  R2 


Figure  5.3:  Roof  and  ridge-line. 


With  an  edge-definition  .  as  the  one  of  DAVIS  (1473)  one  will 
have  difficulty  identifying  an  edge  in  figure  5.3.  DAVIS 
(1973)  defines: 

"An  edge  is  the  boundary  between  two  regions  R^,  R2  of 
different,  constant  gray  values." 

This  could  be  written  as: 

(p,q)  €  *  p  €  Rj_  &  q  e  R2,  f)  r2=0  (5.18) 

If 


3  a<b<c  €  G • (  V 

x1  £  R^ 


V  a  <g(x- )  <  b  <g(x_)  <  c)  v 

x2€r2 


(  V  V  a  <g  (x-)  <  b  <g  (x.)  <  c) 

xl€  R1  x2£  R2 


(5.19) 


then  DAVIS  (1973)  must  call  all  (p,q)  C  L  "edges". 


One  can  employ  the  definition  of  edges  for  multiple 
images.  Let  B  with  the  neighbourhood: 

r(x, y, k) :  =  | (n , v, k)  /n-x/  +  /v-y/  =  l|  (5.20) 

where  K  denotes  the  spectral,  channel. 

Then  the  weight  w  of  an  edge  e 


((Xj^.yj^.kj^) ,  (x2,y2,k2))  is  - - 

w(  (x^y^kj.)  ,  (x2,y2,k2).)  :  *  n-1  J  £  (g  (x2,y2,  i) -g  (x^y^i) 


-66- 


We  see  that  the  weight  of  the  edge  is  independent  of 
image  channel  indices  k^kj. 


5.2  EDGE  DETECTION  OPERATORS 


5.2.1  General  principles 


A  primary  element  of  an  edge  operator  is  the  number 
of  pixels  used  to  compute  the  value  of  an  edge.  There 
exists  a  great  multitude  of  algorithms  and  numerous  reviews 
(e.g.  DAVIS,  1973;  PRATT,  1978;  LEBERL  and  KROPATSCH,  1978b) . 
Of  these  methods,  4  are  selected  for  detailed  analysis: 
three  are  linear,  one  is  non-linear. 

A  general  formulation  for  an  algorithm  must  include 
the  case  of  multiple  input  images:  Image.  ,  k  *  1,  ...n. 

Edge  operations  are  performed  in  parallel:  a  "parallel"  al¬ 
gorithm  is  defined  to  perform  operations  Oj^,  0_,  •••0 
where  the  outcome  is  not  affected  by  the  sequence  of  the 
operations . 


procedure  EDGE  (Image,  n,  f.  A) 
begin 

for  all  i,  J  £  B  do 
A(i, j)  »  f(i,j,n, Image) 

end. 


Ficure  5.4:  PDL-description  of  an  edge  detection  procedure. 


Procedure  FDGE  addresses  image  points  b  £  B  and 
surroundings  U*  of  b  of  the  order  k: 

l^lb):  *  jxCB  |  x  CrS(b),  s  4kj 

The  operations  of  EDGE  can  therefore  be  called  "local".  As 
a  result  one  obtains  a  new  (gray-)  value  g. (b)  for  image 
point  b:  1 

gx(b) :  =  0fa  (g(Uk(b))  . 

The  parallel  algorithm  consists  of  performing  opera  - 
tian  0b  on  all  image  points  b.  It  is  linear  if: 

V  V  0(ax+by)  =  a-O(x)  +  b-O(y). 

x,y£g(B)  a,  b  £  R 
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5.2.2  Linear  algor i  thins 

Typical  linear  parallel  edge  detection  operations 
are:  high-pass  filtering,  oriented  numerical  differenta  - 
ticn ,  gradients,  Laplace-operations .  These  methods  can  he 
described  as  an  operation  g  with  a  mask  m,  that  is  a 
(local)  function 

mx:  U(x)  -*■  R 

on  a  surrounding  U(x)  of  an  image  point  x  C  B.  The  value 
m  (b)  of  b  £  U(x)  is  only  defined  by  the  relative  position 
or  b  with  respect  to  x. 

Operation  gm  then  followsas: 


g  (b)  :  =  S  9  (x)  .  m.  (x). 
m  x  £  U(b) 


Sample  masteas  used  within  an  8-neighbourhood  rg  of  a  2- 
dimensional  image  area  are : 


(a) 

-1 

0 

1 

horizontal  forward  difference  on  U1 (x) 

-1 

0 

1 

-1 

0 

1 

(b) 

-1 

-1 

-1 

Laplace-operator  on  U* (x) 

-1 

8 

-1 

-1 

-1 

-1 

(c) 

-1 

2 

-1 

double  horizontal  forward  difference 

-1 

2 

-1 

on  U* (x) . 

-1 

2 

-1 

(d) 

An 

iterative 

use  of  a  type  of  Laplace-operator 

is 

implied  in 

MONTOTO's  (1977)  method: 

mb 

(b)  = 

(card 

(U (b) ) -1) /card  (u(b)) 

m^x)  = 

-1  / 

card  (U(b))  if  x  /  b. 

5.2.3  Non-linear  algorithms 


Numerous  authors  have  proposed  non-linear  edge  detection 
methods  (compare  PRATT,  1978) .  These  are  usually  more  time- 
consuming  than  linear  ones  and  therefore  mostly  limited  to 
UMx)  of  ro. 

O 


/  > 


5.2.4  Methods  investigated 


The  following  four  methods  were  investigated: 

(0P1)  horizontal  difference 

f:  yij  =  ^  xij  “  xi*lj  ^ 

(OP 2)  Laplace  operator 

f:  yij  =  ^  4xij"xi-lj"xij-l”xi+lj‘xij+l/ 

This  corresponds  to  the  mask  in  rQ: 

0-10 
-1  4  -1 

0-10 


(0P3)  operator  of  Rosenfeld  (see  DAVIS,  1973) 


n  .  . 

n 

f  J 

y  =  - 
1  n 

^  (xij-k“xij+k)j 

1 

This  operator  has  been  implemented  with  n  *  3. Pig. 5. 5 
illustrates  the  image  pixels  entering  into  the  operation 
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(0P4)  operator  of  Roberts  (PRATT,  1978) 

T  2  2l1/2 

f:  yij  “  [(Xij  “  Xi+1  j+l}  +  (xi+l,j  "  xij+l)  J 

This  is  non-linear. 


5.3  TEST  DATA  FOR  EVALUATION  OF  EDGE  OPERATORS 


The  six  synthetic  images  of  figure  5.6  were  generated 
to  test  edge  operators  to  define  an  optimum.  The  images  all 
show  simple  geometric  forms  except  for  figure  5.6  (f)  that 
is  a  segment  of  an  actual  LANDSAT  image  (lake) . 


(e)  Curve 


(c)  Diagonal 


Figure  5.6:  Six  synthetic  test  images 
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A  theoretical  edge  reference  image  is  defined  assigning 
each  image  point  y  £  b  a  value: 


with 


gtYiJ  =  max  (g  (x.  . )  -g  (av) ) 
J  k=l , 8  * 


ak  £  r8  ^xij* 


and  x. .  are  the  pixels  from  the  synthetic  test  image,  gray 
values J are  100  and  110.  An  example  of  an  original  image 
section  and  its  edge  reference  image  is  shown  in  fig. 5. 7  (a), (b) 
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Figure  5.7:  Section  of  synthetic  image,  bright  areas  with 
gray  values  lo,  dark  areas  with  0. 

(a)  image;  (b)  edge  image. 


Similary  edge  reference  images  for  figure  5.6  are  shown 
in  figure  5.8. 


The  test  images  were  provided  with  noise: 


noisy  image  *  synthetic  image  *  noise  *  f, 

where  noise  has  mean  0  and  standard  deviation  1;  f  is  a  noise 
factor.  Noise  added  to  images  is  shown  in  figure  5.9. 


5.4  APPLICATION  OF  THE  EDGE  OPERATORS  TO  THE  SYNTHETIC  IMAGE 


All  four  edge  operators  were  applied  to  all  synthetic 
images  with  noise  factors  0, 1, 2, 3, 4. Then  the  result  was 
evaluated  as  follows: 

(a)  thresholding  of  edge  images  creation  of  binary  edge 

images ; 

thresholds  used  were  60,80,100,120,140; 


(a) Horizontale 


(b)  Verticale 


(c)  Diagonale 


(d)  Sloping  line  (e)  Curve 


If)  Kochel-lake 


Figure  5.8:  Edge  reference  images  for  the  synthetic 
data  of  figure  5.6. 


(b)  comparing  position  and  number  of  edges  in  binary 
image  with  those  in  the  reference  image. 

From  this  one  has  R  correctly  identified  edges,  N  non- 
identified  but  existing  edges  with  value  ^  0,  and  F  incor¬ 
rectly  identified  edges  =  0.  Figure  5.10  illustrates  the 
fact  that  the  sets  R,  N,  F  are  not  disjunct.There  are  totally 
K  edges  that  exist,  NK  non-edges  and  B  image  points. 


We  see  that 


card 

(K) 

=  card  (N) 

+ 

card 

(R) 

card 

(NK) 

=  card  (RNK) 

+ 

card 

(F) 

card 

(B) 

=  card  (RNK) 

+ 

card 

(N)  +  card  (R) 

+  card  (F) . 

From  this 

we 

compute 

PR  = 

100  *  cayd  J,Rj 
100  card  (K) 

• 

» 

PRNK 

-  «»*5& 

(RNK) 

(NK) 

and  finally: 

PR  »  card  (NK)  +  PRNK  ±  card  (K) 
-  “  card  (NK)  +  card  (K) 


Figure  5.11  shows  the  M-values  for  one  test-image  (f)  and 
one  operator  (ROBERTS) . 
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Figure  5.11:  Means  M  of  performance  of  edge  operator  of 
Roberts  used  on  synthetic  image  Kochel. 
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Figure  5.12  shows  an  overview  of  all  operators  and  all  test- 
images,  using  the  best  threshold.  Fig.  5.13  compares  clean 
an  noisy  images. 
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Figure  5.12:  Maxima  of  mean  M  of  performance  factor. 
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Figure  5.13:  Maxima  of  mean  M  in  clean  and  noisy  image 

(in  brackets),  using  maximum  noise  factor  of  4. 


From  all  the?e  results  we  can  conclude 


-  thresholds  have  great  importance; 

-  operators  of  Rosenfeld  and  Robers  are  superior; 

-  the  simpler  operators  depend  strongly  on  the 
direction  of  structures  in  the  image. 


We  also  found  that  with  increasing  noise  the  simpler 
algorithms  improve:  this  is  caused  by  giving  correctly 
identified  edges  a  higher  weight  in  equation  (5.22)  than 
wrongly  identified  ones.  With  noise  more  correct  edges 
are  found,  but  also  more  incorrect  ones.  The  better  edge 
operators  identify  all  edges  with  and  without  noise,  but 
with  noise  they  also  produce  wrong  edges. 


Computing  times  for  the  4  investigated  edge  opera¬ 
tors  compare  as  follows: 


Rosenfeld  Roberts 


Laplace 


horizontal 

difference 


3.75  2.27  2.15  1 


Since  ROSENFELD  and  ROBERTS  are  practically  equivalent 
in  their  overall  performance  it  can  be  concluded  that  the 
method  of  ROBERTS  is  the  most  appropriate  due  to  more 
modest  computing  requirements. 


5.5  OTHER  PRE-PROCESSING  FUNTIONS 


Numerous  other  techniques  must  be  used  for  pre-processing 
prior  to  pattern  analysis.  The  standard  one  are  being  used  such 
as 


-  contrast  manipulation; 

-  principal  component  transformation; 

-  ratioing  of  2  images ; 

-  generation  of  vector  length  images; 

-  clustering  technique  using  principal  component  transform. 

In  addition  a  specific  technique  was  developed  to  enhance 
features  (areas)  in  a  specific  section  of  the  image. 

From  external  sources  the  approximate  size  and  location 
of  a  feature  is  known.  In  the  identified  image  area  a  histo¬ 
gram  is  produced  in  each  image  channel.  Since  the  information 
on  the  feature  is  approximate  only,  the  histogram  will  contain 
pixels  not  belonging  to  the  feature. 

To  suppress  the  effect  of  those  pixels  that  do  not  belong 
to  the  feature,  one  now  selects  a  portion  [A,B]  of  the  histo¬ 
gram  acc.  to  figure  5.14  so  that  a  certain  percentage  of  the 
pixels  is  within  [a,B]. 
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Figure  5.14:  Range  [  A,  B]  contains  80  X  of  all  values  of 
the  histogram. 


Let  us  now  denote  with  UGv»  k  ■=  1,  ..  n  the  lower 
boundary  of  gray  values  (a)  and  OGfc,  k  *  1,  ...  n  the 
upper  one.  k  is  the  index  for  the  channel. 


We  set 


n 

MINP  =  n  UGV  (5 

k=l  * 


n 

MAXP  =  n  OG.  (5 

k=l  * 


and  obtain  the  transformation 


MINP  * 


-) 


_ 255 _ 

MAXP-MINP 


(5 


for  all  image  points  B 


ij 


(k) 


.23) 


.24) 


.25) 


This  transformation  causes  the  image  points  with  gray 
values  similar  to  those  of  the  feature  used  for  computing 
MINP,  MAXP  to  fall  in  the  gray  value  range  (0,255).  All  other 
points  will  be  outside  that  range  (either  0  or  255) .  This 
procedure  works  well  with  large  areas  such  as  lakes, forests, 
etc. .  With  linear  features  (roads)  one  has  to  transform  with 
another,  large  object  and  then  consider  the  0  or  255  gray 
values . 

Figure  5.15  shows  4  original  LANDSAT-channels,  5.16  a 
vector  length  image,  5.17  the  two  first  principal  components 
PCI,  PC2  and  5.18  a  compressed  image  acc.  to  equation  (5.25) . 
For  this  river  Mur  was  used  to  compute  equations  (5.23), 

(5.24)  and  this  causes  the  free-way  to  fall  into  gray  value 
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Figure  5.16:  Vector  length,  image  of  scene  in  figure  5.15 
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6.  RECOGNITON  PROCEDURES 


Recognition  in  our  special  case  is  concerned  only  with 
separate  objects.  These  are  objects  which  appear  in  a  map 
and  have  either  a  uniform  greyscale  or  a  uniform  texture. 
Images  with  objects  of  uni form  texture  can  be  transformed 
into  images  of uniform  greyscale  by  one  of  the  methods  e.g. 
described  by  WESZKA  (.1979.)  where  the  resulting  grey  values 
are  texture  measures  of  a  local  environment.  So  we  can  con¬ 
centrate  on  features  in  images  of  uniform  grey-scale. 

A  second  restriction  to  general  recognition  procedures 
comes  from  the  fact  that  we  have  already  a  lot  of  informati.n 
from  the  map  concerning  the  object:  we  know  the  approximate 
location  of-  the  object,  we  know  its  shape,  possible  the  range 
of  grey  values.  It  is  therefore  meaningful  to  start  with  a 
hinary  matrix  of  the  shape  of  the  object  to  use  it  as  a  mask 
to  compare  the  appropriate  grey  values  with  the  predicted 
value  from  the  data  bank  or  from  other  calculations  on  the 
image.  The  result  of  the  comparison  is  then  another  binary 
matrix  called  "  pattern  ".  These  parameters  of  the  first 
step  of  the  recognition  process  are  also  shown  in- figure  2.7. 

In  a  second  step  one  has  to  decide  whether  the  pattern 
found  in  step  1  represents  the  object  or  not.  Here  it  is 
called  VERIFICATION.  It  computes  a  measure  of  similarity  bet¬ 
ween  object  and  pattern.  An  object  is  denoted  as  rejected 
(not  recognized)  when  the  verification  measure  is  below  some 
accept  treshold  value. 

Four  methods  for  centering  the  map  feature  correctly 
over  the  image  have  heen  studied  (see  figure  2.7).  Their  prin¬ 
ciples  of  operation  and  implementations  are  described  in  the 
following  chapters,  together  with  the  verification  functions 
employed.  These  chapters  are  preceded  by  definitions  of  cor¬ 
relation.  functions,  which  are  used  in  several  of  the  proce¬ 
dures.  Corresponding  practical  experiences,  results  and  tests 
are  presented  in  a  separate  chapter  8  so  that  the  current 
chapter  is  merely  theoretical. 


6.1  CORRELATION  FUNCTIONS  WITH  BINARY  MATRICES 


From  the  registration  of  an  image  with  another  image  we 
know  many  correlation  functions  (DOWMAN,  1977;,  BARNEA  et  al., 
1972;  ANUTA,  1970) .  Correlation  functions  with  binary  matrices 
are  treated  here  because  the  formulas  restrict  themselves  to 
expressions  which  use  only  standard  statistical  measures  and 
masked  image  areas.  This  advantage  has  two  effects:  the  number 
of  accesses  to  the  image  elements  can  be  minimized  and  standard 
procedures  for  statistical  measures  can  be  used. 

Our  formulas  are  derived  from  two  different  functions : 

The  coefficient  of  correlation  (KREYSZIG,  1973)  and  the  cross¬ 
correlation  (BARNEA  et  al.,  1972). 


A 
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The  following  statistical  standard  measures  are  used  in 
our  formulas  (binary  matrix  and  image  S  should  have  the  same 
window  size) : 


(6.1)  p (m)  =A  M(i,  j)  . surface  of  binary  matrix  M 

i/ j 

(6.2)  F(S)  =  n  *  m  .  surface  of  an  image  window  S  of  n 

columns  and  m  rows 

(6.3)  ^.(S)  =  (XT  S  (i,  j)  )/F{S)  . . . .  mean  value  of  image  window  S 

^  i*  j 

(6.4)  ^(M)  =  F  (M)/F(S)  ...  mean  value  of  binary  matrix 


(6.5)  *cM(S)=[^(S(i<j)*M(i<  j))J  /F 

binary 

2 

(6.6)  Q(  S)  =[17(3(1,  j)-  MS))2]  /(F(S)-l) 

i » j.  ^ 


(M)  ... 
binary  mask  M 


mean  of  image  S  over  the 

. . .  variance  of  image  S 


6.1.1  Derivation  of  the  coefficient  of  correlation 
n  m 

2T  ZZ  (S(i'j)  -£<S))*<M(i.J>  -/U/(M)) 

_  .1=1  fel _ 


coco 


]/  [±  g«U.J)  -  £•«> )  2J  *[£ 


Building  the  product  and  splitting  the  sums  simplifies  the 
numerator  rN: 


n  m 


n  m 


n  m 


rN= 


ZZZZ  S(i,j)M(i,j)-|U(S)ZZ^  M(i,  j)-£c(M)]^  ZH  S(i,j)+^as^(M)F(S) 


Using  equations  (6. 1) ,  (6. 3) ,  (6 . 4)  and  (6.5)  from  above  rN  becomes 
rN=  F(H).f[?(S)-f(S)F(H)-F(M)ms)*  F(S)  £.(S)+£.(S)|j|j-  F(S)  . 

It  follows  directly 

rN=  F(M)£/M(S)-F(M)  £t(S)-F(M)  £t(S)  +F  (M)  £  (S) 

F(M)^m(S)- £-(S)). 


and 


rN= 
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1 


With  equ.  (6.6)  and  knowing  that  M2 (i, j)=M(i, j)  the  square  rD^ 
of  the  denominator  rO  becomes 


rD2=  (F(S)-l)*  ^2(S)*^ZT  M  ( i  /  j )  -  F(S)  (M)j 

and  from  (6.1)  and  (6.4)  it  derives 

rD2=  [f (s)  -lj  *  er2(s)*  [f(m)-  f(s)*f2(m)/f2(s)] 

and 

rD2*  [F(S)-l]  *  S^(S)*  F(M)/F  (S)*[f(S)-F(M)J  . 


Hence  the  quotient  of  rN  and  rD  results  to 
,  - 


(6.7)  coco 


fr(S) 

laTsj 


[  F(S) 


F(M)  F( 
)-  lJ  *fF 


IS) _ 

(S)-F(M)J 


In  an  analogous  way  the  formula  for  the  cross-correlation  derives 
from 


cross= 


7^  M(i,j)*S(i, 


j) 


]/[£  M2(i.j)j»[2i  s2u.j>j 


and  results  in 

(6.8) 


cross  =  £*m(S)  * 


F  (M) 


F(S)  ^-2  (S)  -l-fF  (S)  -l]  0’2(S) 


6,1,2  Correlation  between  two  binary  matrices  M, ,  M, 

-  •  '  '  .  4 


In  the  case  of  two  binary  matrices  M,  and  M_  only  surfaces 
F(M)  and  F(S)  are  of  interest.  M  is  one  or  the  two  matrices  or 
their  product  M-.M,,  and  F(S)  represents  the  minimal  window 
size  which  comprises  M.  and  M2*  Pixels  outside  of  the  window 
of  one  of  the  binary  matrices'are  per  definition  set  to  0. 

Since  the  derivation  of  the  formulas  is  trivial  only  the  re¬ 
sults  are  presented  here. 
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F(S)*F(M1.M2)  -  F(M1)#F(M2) 


(6.9) 


bicoco  = 


Y  FfM^aFfM^ljtS)-  F(Mx)j  *£f(S)-  F(M2)]' 
Binary  cross-correlation: 


(6.10) 


F(M1  .  M2) 


bicross= 


y  F(M1)^  F(M2) 


For  comparison  purposes  formulas  (6.7)  and  (6.9)  of  the 
coefficient  of  correlation  are  used  as  absolute  values.  All 
correlation  measures  which  appear  in  the  following  chapters 
are-  given  in  percent.  That  means  that  the  values  are  normalized 
between  0  and  100,  where  0  stands  for  no  correlation  and  100 
for  the  best. 


6.2  COMPARING  THE  REAL  WITH  A  SYNTHETIC  IMAGE 


It  is  obvious  that  the  binary  image  constructed  from  the 
map  feature  is  a  mask  or  a  template  or  a  synthetic  image.  This 
could  be  shifted  over  the  window  to  find  maximum  corresoondence 
between  the  window  pixels  and  the  template.  This  follows  the 
sequential  similarity  detection  algorithm  of  BARNEA  and  SILVERMAN 
(1972)  which  is  normally  applied  to  correlating  two  images. 
However,  the  method  is  typically  suited  to  define  a  geometric 
distortion  vector;  it  is  not  equally  suitable  to  define  the  shape 
of  a  feature  in  the  image. 

One  of  the  possible  implementations  of  this  shift  algorithm 
is  shown  in  figure  6.1. 

It  starts  at  the  position  of  the  binary  matrix,  which  is 
estimated  from  the  map.  Correlation  measures  max0  are  calculated 
for  levels  of  equal  distance  d  from  the  original  positicr  It 
continues  for  ascending  distances  d  until  either  no  better 
correlation  is  found  or  a  maximum  distance  Dmax  is  reached.  The 
function  binary  matrix  (pos)  delivers  the  binary  matrix  shifted 
to  position  "pos".  The  condition  in  line  10  of  figure  6.1  defines 
all  positions  of  fixed  distance  d  from  the  original  "position" 
in  the  8-neighbourhood  (figure  6.2).  The  position  "positmax"  of 
the  best  match  is  used  to  shift  the  binary  matrix  to  the  resulting 
pattern . 


(1)  procedure  SHIFT  (image,  binary  matrix,  pattern) ; 

(2)  begin 

(3)  max0  : =0 ; 

(4)  maxi  :=  100  *  Abs (coco  (image  ,  binary  matrix)); 

(5)  d: ®  1; 

(6)  positmax:=  position (binary  matrix); 

(7)  while  maxl>  max0  or  d^Dmax  do 

(8)  begin 

(9)  max0;=  maxi; 

(10)  for  all  pos  with  (max (x (pos-position) , y (pos-position) ) =d) do 

(11)  begin 

(12)  pattern :=  binary  to  matrix  (pos) 

(13)  c ;  =  100  *  Abs (coco  (image, pattern) ) ; 

(14)  if  c>  maxi  then 

(15)  begin 

(16)  maxl:=c; 

(17)  positmax:=  pos; 

(18)  end; 

(19)  end; 

(20)  d:=  d+1 ; 

(21)  end; 

(22)  pattern:=  binary  matrix (pos itmax) ; 

(23)  end. 


Figure  6.1:  Recognition  procedure  SHIFT. 
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Figure  6.2:  Distances  in  the  8 -neighbourhood 


6.3  REGION  DETECTION  BY  THRESHOLDING 


The  principle  of  the  SHIFT-method  was  to  locate 
the  real  position  of  the  object  in  the  image  by  corre¬ 
lating  the  binary  mask  and  the  image  at  different  hori¬ 
zontal  and  vertical  positions.  The  shape  of  the  object 
never  changed,  so  that  local  distortions  could  not  be  de¬ 
tected.  We  tested  another  method  of  region  recognition, 
which  starts  with  a  rudimentary  binary  feature  and  assemj 
les  pixels  belonging  to  a  region  R  according  to  the 
corresponding  grey  values  in  the  picture. 

The  simplest  method  for  assembling  the  pixels  x.  be¬ 
longing  to  a  region  R  is  by  thresholding.  Supposing  Chat 
we  start  with  the  mask  M  from  the  map  data  bank  projected 
on  the  image  we  calculate  the  following  parameters  before 
thresholding  the  image: 

-  a  lower  and  an  upper  threshold  value  and 

-  a  starting  mask  S  containing  all  pixels  which 
are  possible  candiates  to  belong  to  region  R. 


If  we  assume  a  dislocation  of  any  pixel  with  a  maximum  of  0 
pixel  diameters  the  starting  mask  S  can  be  computed  from  M 
by  adding  a  band  of  D  pixels  around  the  feature.  Then  for  all 
pixels  of  R,  S(i,j)  =  1  under  that  assumption  (see  fig. 6. 3b). 

T'o  :  choose  the  threshold  values  we  tested  various  me¬ 
thods  with  different  complexity.  The  simplest  and  fastest 
one  is  to  take  prefixed  values.  That  is  possible  if  the  un¬ 
derlying  image  has  been  preprocessed  so  that  all.pixels  be¬ 
longing  to  a  certain  object  type  have  a  known  range,  which 
is  then  used  as  threshold  interval  £lPU,  IPO]  . 

Sometimes  we  know  only  that  the  pixels  of  the  object 
accumulate  at  one  end  of  the  grey  scale.  In  this  case  a  hi¬ 
stogram  can  be  used  to  find  the  other  value.  The  histogram 
shall  be  taken  from  an  image  region  which  contains  the  whole 
object.  This  could  be  a  rectangular  region  (S-window  in 
figure  6.3a)  D  pixels  greater  in  each  axis  direction  than  the 
original  window  of the  object  or  the  starting  mask  S  as 
described  above.  Since  the  real  surface  of  the  object  is  known 
from  our  data  bank  the  threshold  value  is  defined  by  the  condi¬ 
tion  that  the  number  of  grey  values  between  the  threshold  pi¬ 
xels  equals  the  surface  content  of  the  object.  In  this  way 
procedure  ADAHIS  works  in  figure  6.3a.  Factor  p  depends 
upon  the  dislocation  0  and  the  shape  of  the  object.  It  is  di¬ 
rectly  proportional  to  D  and  to  the  ratio  F(S)/F(M)  where  M 
stands  for  the  binary  matrix  representing  the  object  and  S 
for  the  starting  mask.  Ftx)  denotes  the  surface  content  of 
the  respective  parameter  x. 

If  no  information  of  the  grey  values  of  an  object  can 
be  predicted,  standard  statistical  measures  produce  sometimes 
acceptable  results.  We  calculate  the  mean  value  MV 
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(1)  procedure  THRESH  (image,  M,S) 

(2)  begin 

(3)  S-window:*  Window  (M)  *  (-D, -D, -fc2D, -fe2D)  ; 

(4)  STARTINGMASK  (binraat , S-Window,  D,S) ; 

(5)  HTSTO  (image,  S,  S-Window) ? 

(6)  ADAHIS  (peP(M),  IPU,  IPO); 

(7)  for  all  (i, j) £  S-Window  do 

(8)  if  S  ( i ,  j )  =  l&-t(IPU  ^  image (i , j ) ^IPO)  then  S(i,j):«fib 

(9)  end. 


Figure  6.3a:  Recognition  procedure  THRESH. 
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procedure  STARTINGMASK Cbinmat ,  window, D, mask) 
begin 

for  all  (i,j)£  window  do 
begin 

if  binmat (i , j ) =1  then  mask  (i,j):=  1  else 
begin 
d:-l; 

sum: =0 

while  sura=0  &  d-D  do 
begin 
k:=d; 

1:=0; 

while  k  0  &  sum  =  0  do 
begin 

sum:  =sum-fchinmat  (i+k,  k-tl)  + 
binmat  (i-1 ,  j  -fek)  + 
binmat (i-k, j-1)  -fc 
binmat  (i+l,j-k); 

k : =k- 1 ; 

1:=1+1; 

end 

d: =d+l ; 
end; 

if  sum=0  then  mask(i,j):*0  else  mask(i, j) :®1; 

end 

end 

end. 


Figure  6.3  b:  Procedure  STARTINGMASK 


and  the  standard  deviation  SO  of  the  grey  values  in  a  certain 
region.  Then  the  threshold  values  are  IPU:»MV-SD  and 
IPO:»MV+SD  respectively.  The  underlying  region  can  be  the  project¬ 
ed  object  matrix  M  or  the  starting  mask  S  or  a  region  smaller 
than  M  or  an  image  window  containing  M  or  only  a  portion  o  M. 

After  these  preparations  the  proper  thresholding  is  per¬ 
formed.  Pixels  from  S  are  eliminated  if  the  corresponding  grey 

/ 
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value.  does  not  fall  into  the  threshold  interval  (.see  lines 
7  and  8  of  figure  6.3a}. 

The  disadvantage  of  the  method  THRESH  is,  that  the  re¬ 
sulting  hinary  mask  S  is  not  necessarily  a  connected  set.  An 
advantage  is  the  parallel  operation  of  the  algorithm. 


6.4  ADAPTING  THE  HISTOGRAM 


Another  method  which  takes  into  account  the  disadvantage 
of  thresholding  starts  with  the  contour  of  the  object  as  gi¬ 
ven  in  the  map.  A  histogram  of  the  image  over  the  binary 
region  will  reveal  that  pixels  are  included  that  should  be 
eliminated.  On  the  other  hand,  pixels  in  the  remainder  of  the 
image  are  not  part  of  the  approximate  binary  projection,  of 
the  object,  but  are  adjacent  to  it  and  should  be  added,  because 
their  grey  values  are  very  close  to  that  of  the  feature.  A 
search  along  the  boundary  of  the  feature  will  thus  permit 
to  connect  the  original  contour,  using  the  histogram  to  decide 
which  pixels  to  reject  and  which  ones  to  accept.  This  process 
can  be  repeated  iteratively  until  either  no  corrections  to  the 
contour  are  done  in  the  last  step  or  a  maximum  iteration  count 
is  reached. 

Figure  6.4  shows  the  procedure  in  more  detail.  MH  is 
there  a  binary  matrix  which  guarantees  that  changing  the  value 


(1)  procedure  ADAPT  (image,  M,  binadapt,  maxiter) 

(2)  begin 

(3)  S- Window :=  Window (M)  + {-maxiter, -maxiter,  +2  maxiter,  +2  maxiter) 

(4)  iter:=  1; 

(5)  alt:=-l ;  neu:=  0/ 

(6)  binadapt :=  M; 

(7)  while  iter  <  maxiter  &  neu>alt  do 

(8)  begin 

(9)  alt:=  neu; 

(10)  MH:  =  (1)  -  binadapt; 

(11)  HISTO  (image,  binadapt,  S-Window) ; 

(12)  ADAHIS  (p«F(M),  IPU,  IPO); 

(13)  for  all  pixel  e S-Window  do 

(14)  if  MH  fcixel)  =0  v  (pixel  €  border  (MH)  &  ( IFUslmage  (pixel )  s  IPO)) 

then  binadapt  (pixel) :=  0 

(15)  else  binadapt  (pixel)  :=  1; 

(16)  MH:=(D-  binadapt; 

(17)  for  all  pixel  € S-Window  do 

(18)  if  MH  (pixel)  ^  v (pixel  £  border  (MH)&  not  (IPU^image  (pixel )^IP0)) 

then  binadapt  (pixe]):=  0 

(19)  else  binadapt  (pixel) :=  1; 

(20)  iter:=  iter+1; 

(21)  neu:=  F(|binadapt  -  M|  ); 

(22)  end 

(23)  end 


Figure  6.4:  Recognition  procedure  ADAPT. 
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o f  a  pixel  from  1  to  0  does  not  influence  the  property 
"pixel C border  (MH) "  during  one  iteration  step.  Lines  (10) 
and  (16)  assign  the  binary  complement  of  binadapt  to  MH. 

The  selection  of  threshold  values  IPU  and  IPO  in  lines  (11) 
and  (12)  works  analogous  to  procedure  THRESH.  Condition 
"neu>alt"  in  line  (7)  and  the  calculation  of  "neu"  in 
line  (21)  to  be  the  surface  content  of  the  absolute  binary 
difference  between  the  original  M  and  the  adapted  matrix 
"binadapt"  express  the  termination  of  the  loop  if  no  more 
changes  have  accurred  and  "neu"  is  empty. 

An  advantage  of  this  procedure  lies  in  the  fact  that 
all  changes  of  the  binary  matrix  preserve  its  connectivity. 
The  idea  for  this  method  derives  from  a  skeletonizing  or 
shrinking  algorithm  of  ROSENPELD  (1970).  "ADAPT"  combines 
the  connectivity-preserving  deletion  process  with  an  addi¬ 
tional  threshold  condition  on  the  image.  ROSENPELD  did 
show  that  the  binary  property  of  his  algorithm  preserves 
connectivity.  Since  we  use  only  the  deletable  pixels  from 
ROSENFELD  his  result  must  be  a  subset  of  ours.  Hence  our 
binary  matrix  has  the  same  property. 

A  disadvantage  of  ADAPT  is  its  relatively  high  processing 
time  compared  with  thresholding.  The  algorithm  has  to  inspect 
for  each  pixel  a  3  x  3  pixel  neighborship  of  the  binary 
matrix  to  test  for  a  border  point  (lines  (14)  and  (18))  and 
in  some  cases  the  grey  value  in  the  image.  This  is  done 
twice  (for  adding  pixels  and  deleting  pixels)  with  an  addit¬ 
ional  histogram  building  at  each  iteration.  The  histogramming 
is  not  absolutely  necessary  at  each  iteration  but  it  responds 
better  to  the  fact  that  before  adapting,  a  great  amount  of 
erroneous  pixels  affect  the  histogram  and  the  choice  of  the 
threshold  values  IPU  and  IPO.  In  the  version  of  figure  6.4 
these  are  updated  at  each  step  in  order  to  get  the  most 
homogeneous  region  in  the  image. 

From  these  considerations  follows  that  ADAPT  should  be 
used  for  the  recognition  of  objects  with  small  deviations 
from  the  expected  location  in  the  image  to  get  a  registra¬ 
tion  of  greater  precision  than  with  other  methods  like 
threshold. 


6.5  LINE  DETECTION 


Lines  in  satellite  images  appear  either  as  borderlines 
between  two  distinct  regions  or  as  linear  features  like 
rivers,  roads  etc.  Since  borderlines  are  produced  implicitly 
when  their  adjacent  regions  are  recognized,  their  recognition 
is  not  of  special  interest  for  us .  On  the  other  hand  they 
do  not  appear  in  our  map  data  bank  as  separated  objects.  However, 
linear  features  do  appear  as  skeletons  in  the  map  because 
they  are  considered  to  be  very  small  regions . 

Several  procedures  which  detect  lines  in  digital  images 
are  proposed  in  the  literature.  In  principle  there  are  3  ways 
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to  find  lines  from  which  the  course  is  known  approximately: 

1)  Shifting  and  rotating  a  binary  mask  over  the  image 
and  determining  the  correct  position  by  the  maximum 
of.  a  correlation  function  like  method  SHIFT  in  chapter 
6.2.  HACK  (1975)  and  SWEDLOW  (1978)  correlated  two 
binary  edges  images. 

2)  Thinning  a  binary  region  produced  from  the  map  projec 
tLon  of  the  skeleton  followed  by  a  thickening  like  the 
production  of  the  starting  mask  in  method  THRESH.  A 
thinning  algorithm  was  proposed  by  ROSENFELD  and  DAVIS 
(1975).  The  thinning  algorithm  eliminates  i n  our  case 
the  pixels  according  to  their  grey  values  like  pro¬ 
cedure  ADAPT. 

3)  Following  the  lines  sequentially. 

The  first  two  possibilities  can  be  performed  by  procedures 
THRESH  and  ADAPT.  We  therefore  concentrate  in  the  following 
on  sequential  line  following. 

Literature  (VANDERBRUG,  1973;  DAVIS,  1973;  HOLDERMANN,  1971; 
TENENBAUM,  1978)  of  sequential  line  following  algorithms 
shows  a  great  similarity  with  general  problem  solving  algo¬ 
rithms  in  Artificial  Intelligence  (NILSSON,  1971;  VEILLON , 

1974).  The  algorithm  in  figure  6.5  represents  a  general  line 
following  algorithm  as  a  result  of  adapting  the  general  pro¬ 
blem  solving  algorithm  of  VEILLON  (1974)  to  the  special  pro¬ 
blem  of  finding  lines.  All  sequential  line  following  algo¬ 
rithms  that  we  found  in  the  literature  are  expressable  by 
this  general  procedure  simply  be  selection  of  special  para¬ 
meters  . 

The  algorithm  LISU  in  figure  6.5  works  on  a  general  search 
space  X.  In  our  case  this  would  be  either  a  window  of  pixels 
containing  the  line  or  the  set  of  pixels  defined  by  the  start¬ 
ing  mask  S  as  used  in  procedure  THRESH  (chapter  6.3) .  The  search 
starts  at  pixel  x  in  the  search  space  X.  This  pixel  must  be 
determined  approximately  by  the  map  projection  of  the  line. 
Important  parameters  are  the  cost  function  c(x,y)  defined  on 
pairs  x,y  of  neighbour  pixels  and  the  heuristic  function  h'  (x) . 
They  essentially  direct  the  search  to  a  aoal  pixel  g£  X  for 
which  "end  (g) "  is  true.  The  heuristic  h' (x)  estimates  the 
sum  of  costs  for  the  optimal  path  from  the  argument  pixel  x  to 
the  goal. 

The  state  of  the  search  is  described  by  a  set  S  of  pixels 
from  X  and  the  approprate  values  ^  and  pointers  "last".  After 
each  step,  (^ (x)  contains  the  minimal  sum  of  costs  for  the 
optimal  patn  from  x  to  x.  The  path  (x  , x,...x.  ,  x)  is  defined 
recursively  by  "las®": 


last  (x)  *  xk  and  last  (x^  =  xi_i  for  i  *  1#  ...,k 

and  the  set  S  contains  all  endpoints  of  minimal  paths  or  in 
other  words,  all  candidate  pixels  to  continue  to  find  the  end 


(1) 

procedure  LISU  (xQ,X,  end,  c. 

h',xn. 

(2) 

begin 

(3) 

xn5=  xo  * 

(4) 

last  (xQ) :=  nil; 

(5) 

for  all  x  £  X  do  6>(x):  =  o®; 

(6J 

( p(xQ ) :=  0; 

(7) 

S:  =  {  x-j  ; 

(8) 

while  t  end  (xn)  4  S  j*  0 

do 

(9) 

begin 

(10) 

eliminate  xn  from  S; 

(ID 

for  all  x  €  (  P  <x_)  -  P  (last  (xn) 

(12) 

begin 

(13) 

(x„)  +  c  (xn , 

x)  ; 

(14) 

if  tpn  <<*>(x)  then 

(15) 

begin 

(16) 

9(x)  :=  <Pn; 

(17) 

enter  x  in  S; 

(18) 

last  (x)  :=  x„ 

(19) 

end 

(20) 

end 

(21) 

m:=  co  ; 

(22) 

for  all  x  6  S  do 

(23) 

if  (p(x)  +  h'  (x)  < 

m  then 

(24) 

begin 

(25) 

m:=  Cf?(x)  +  h'  (x)  ; 

(26) 

xn:-  X 

(27) 

end 

(28) 

end 

(29) 

end 

Figure  6.5:  A  general  line  following  algorithm. 


pixel.  These  invariance  conditions  hold  also  for  the  goal  state 
Therefore  the  successfull  result  is  defined  by  the  pointer 
chain  "last",  starting  at  the  goal  pixel  x  and  the  minimal 
cost  <f7(xn)  at  that  pixel.  Condition  S=0  can  be  true  only  in 
the  case  of  an  error  in  the  definition  of  X  or  if  for  no  pixel 
x  of  X  end  (x)  takes  the  value  "true". 

The  importance  of  the  heuristic  lies  in  its  tendency  to 
prefer  paths  which  tend  towards  the  goal.  Since  t0(x)  contains 
the  cost  sum  for  the  path  from  x  to  x  and  h'  lx)  estimates  the 
cost  sum  for  the  rest  of  the  path  from  x  to  a  goal  pixel, 

V7  (x)th'(x)  represents  in  line  (23)  an  estimate  for  the  path 
passing  through  x.  The  pixel  x  of  S  where  the  search  shall 
continue  at  the  next  step  is  selected  in  lines  (22)  until  (27) 
to  be  on  the  path  for  which  minimal  total  costs  are  expected. 


This  informal  description  of  the  algorithm  is  based  on  the 
proof  which  can  be  found  in  VEILLON  (1974) . 
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6.5.1  Cost  and  heuristic  function  for  some  special  methods 


One  method  often  mentioned  in  the  literature  is  the 
"heuristic  search"  of  a  line.  Cost  and  heuristic  functions 
c(x,y),  h' (x)  are  determined  according  to  gray  values  g(x) 
and  g(y) .  MARTELLI  (see  VANDERBRUG,  1973)  uses  as  search 
space  X  the  set  of  all  possible  edges  on  the  image  B.  The 
set  of  neighbors  on  X  is  defined  as  follows: 

If  a,b,c,d,e  and  f  denote  points  of  image  B  in  one 
of  the  following  constellations , 

c  d 

a  b  or 

f  e 

then  the  neighbors  of  edge  (a,b)  are  enumerated  in: 

P((a,b)):=  £(a,c)  ,  (a,f)  ,  (b,d)  ,  (b,e)  ,  (c,d)  ,  (e,f)  J 
The  cost  function  is  calculated 

c ( (a ,  b)  ,  (x,y))  :  —  max  {0,  2M-  (g  (a)  +g  (x)  )  +  (g  (b)  +  g(y) )} 

with  M  to  be  the  contrast  g(p)-g(q)of  an  optimal  edge  (p,q) . 
The  search  starts  with  an  edge  x  =  tll,s  ),  (l,s  *1))  behind 
column  s  in  the  first  line  of  the  picture  and  ends  when 


f  a  c 
e  b  d 


end  ( (a,b) , (x,y) ) :=  (a=last  line) & (x=last  line) 
becomes  true. 

The  algorithm  uses  further  the  following  heuristic: 

r  0  if  line  (a)  >  max  {line  (x)f(x,y)£S^  -3 

h' ( (a, b) ) :  =  J 

\  m  else 

The  heuristic  prevents  the  search  to  follow  paths  from 
end  edges  which  lie  3  lines  before  the  highest  line  number 
which  had  ever  been  developed  and  stored  in  set  S  during  the 
preceding  search.  In  effect  backtracking  is  reduced  and  the 
search  is  forced  towards  the  goal,  the  last  line. 

The  Look-Ahead  algorithms  and  the  algorithm  of  CHERNIAVSKY 
(VANDERBRUG,  1973)  use  a  similar  cost  function  and  estimate 
the  heuristic  h' (x)  in  some  neighbourhood  U(x). 


An  other  efficient  group  of  search  methods  are  summarized 
under  the  notion  "Directed  Search".  The  search  is  there  directed 
by  external  informations,  which  cannot  be  calculated  directly  on 
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the  image  and  come  therefore  from  the  outside.  They  determine 
the  heuristic  h‘ (x)  and  influence  via  h'  the  flow  of  the 
search.  A  special  form  of  the  "Directed  Search"  is  called 
"Planning"  (KELLY,  1971;  DAVIS,  1973) .  It  assumes  that  inter¬ 
mediate  goals  are  known  from  preceding  steps  and  uses  them 
through  the  heuristic  to  follow  the  line  passing  through. 

HARLOW  (in  DAVIS,  1973)  describes  an  algorithm,  which 
gets  its  information  directing  the  search  from  a  description, 
of  the  scene.  This  early  approach  could  be  regarded  very  near 
to  ours  in  which  a  map  knowledge  shall  help  to  find  lines  in 
satellite  images. 

A  method  proposed  by  MONTANARI  (VANDERE.RUG,  1973),  which 
falls  into  the  category  of  "Dynamic  Programming",  is  in 
reality  a  simplification  of  the  general  algorithm  by  setting 
the  heuristic  h'  equal  to  zero  for  all  arguments.  Further 
restrictions  arise  for  the  set  of  neighbors  allowing  only 
neighbors  in  the  original  sense  for  which  certain  conditions 
for  their  grey  values  are  fulfilled.  The  search  space  is  here 
built  by  the  image  pixels  B  itself. 

A  similarity  between  the  "Bidirectional  Search"  (VANDER- 
BRUG,  1973)  and  the  scheme  of  the  general  sequential  line 
following  algorithm  is  not  obvious .  The  search  starts  there 
simultaneously  from  start  and  end  point  of  the  line.  For 
the  starting  point  a  set  S  of  alternatives  is  created  and 
for  the  end  point  a  set  T  for  paths  running  into  the  goal. 

A  difficulty  arises  for  the.  computation  of  the  "end “-function 
which  has  the  form: 

end  (x)  <=>  x€SaT 

This  derives  from  the  fact  that  the  function  must  be  calcu¬ 
lated  at  each  step  and  requires  an  inspection  of  every 
element  of  S  and  T. 

HOLDERMANN  proposed  in  1971  a  straightforward  method. 

His  algorithm  selects  edges  in  the  search  space  giving  them 
two  values:  a  weight  and  a  direction.  The  neighbourhood  p(x) 
has  the  form  of  a  sector  centered  in  x  of  several  pixel 
radius.  Its  main  direction  is  the  same  as  the  direction  of 
edge  x.  The  characteristics  of  the  cost  function  must  be 
selected  so  that 

1)  the  costs  for  a  chosen  alternative  are  always  lower 
than  for  other  alternatives, 

2)  a  selected  edge  corresponds  best  to  its  predecessor 
in  weight  and  direction  and 

3)  if  there  is  no  direct  successor,  adjoining  sectors 
are  selected. 


The  heuristic  has  no  importance  in  this  approach. 


UNCLASSIFIED 
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6.5.2  Interrelationship  between  cost  and  heuristic  functions 


VEILLON  (197^  enumerates  2  conditions  for  costs  c  and 
heuristic  h1: 


(A)  c(x,y)>0  for  all  x,y  of  X  and  xe  T(y) 

‘n-1 


(B)  0  £  h'  (x)  S  min]  ^  c(xi,xi+l) 

i= 


x=X2.  &  end(xn)  &  xi+1€ 


Hx,) ,  i=l ,  n-lj 


In  the  case  of  line  detection  the  costs  c  depend  essen¬ 
tially  on  the  grey  values  of  the  corresponding  pixels.  If 
we  assume  that  the  image  has  been  preprocessed  so  that  pixels 
on  the  line  have  high  grey  values,  then  the  desired  path  has  the 
property  that  the  sum  of  grey  value  differences  to  the  maximum 
is  minimal.  This  leads  to  the  conclusion,  that  c(x,y)  depends 
only  on  the  grey  value  g(y).  We  further  restricted  the  functional 
class  of  c  to  the  form 


c(x,y)  ■  cg(g(y))  =  (a  +b*g(y))n  with  a,b,neR  and  n>0. 


The  choice  of  c  was  a  result  of  extensive  tests  and  the  con¬ 
sideration  that  pixels  near  the  optimal  grey  value  should 
yield  small  costs  while  grey  values  far  away  from  the  optimal 
one  should  be  considered  only  in  the  worst  case.  Parameters-  a, 
b  and  n  are  determined  by  3  conditions  which  cam  be  derived 
from  our  data  base  knowledge  if  we  assume  that  h' (x)  measures 
the  Euclidean  distance  from  the  goal. 


1.  If  B  is  the  highest  grey  value  of  the  search  space  X 
then  we  set 


(Bl) 


eg (B  +1)  =  0. 


2.  Since  pixels  of  the  line  accumulate  at  B,  we  choose  a 
range  [C,B]  of  grey  values  in  which  all  line  pixels 
should  lie.  Let  H(g)  denote  the  histogram  of  the  grey 
values  g  in  the  search  space  X.  Then 
B 


y  H(g)  =  L  can  be  used  to  determine  C  if  L  is  the 


g=c 


length  of  the  line  in  units  of  pixels .  To  guarantee 
general  condition  (B),  we  demand  that  eg  (x)  for  x£^C,Bl 
should  be  1  in  the  mean  and  therefore: 


B 


(g)  *  eg  (g) 


(B2) 


SES. 


=  1 


3.  Pixels  outside  of  the  line  should  get  high  costs.  We 
set,  the  mean  costs  for  such  pixels  x  with  g(x)£[A,c  -  i] 
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equal  to  the  total  sum  L  over  all  goal  pixels. 

"A"  stands  here  for  the  minimal  grey  value.  That  means: 


<B3) 


H  (g)  *  eg  (g) 


iSH(9’ 


L  . 


From  these  3  condidtions  parameters  a,b  and  n  are  derived  de¬ 
livering  a  consistent  pair  of  costs  and  heuristics  (cg,h'). 


6.5.3  Representation  of  the  result  of  Line  following 


The  result  of  procedure  LISU  consists  in  a  pointer  chain 
"last"  from  a  goal  pixel  x  back  to  the  starting  point, .  The 
transformation  in  a  binarymatrix  "binline"  is  done  in  lines 

(8)  to  (13)  of  procedure  LINDET  (figure  6.6) .  Procedure 
COSTPARAMETER  determines  the  values  of  a,  b  and  n  with  the 


(1)  procedure  LINDET  (image,  mapline,  start,  goal,  binline) ; 

(2)  begin 

(3)  S-Window:=  Window  (mapline) + (-D.-D, +2D, +2D) ; 

(4)  STARTINGMASK  (mapline,  s-window,  D,  Snask) 

(5)  HISTO  (image,  Smask,  S-Window) 

(6)  COSTPARAMETER  (a,  b,  n) ; 

(7)  LISU  (start,  Smask,  end,  c,  h' ,  xn,  last) ,• 

(8)  for  all  x  e  S-Window  do  binline  (x)  :=*  0; 

(9)  while  xn  nil  do 

(10)  begin 

(11)  binline  (x_) :=  1; 

(12)  xn:=  last  (xn) 

(13)  end 

(14)  end. 


Figure  6.6:  Filling  a  binary  matrix  "binline"  with  pixels 
of  a  line. 


aid  of  the  histogram  built  in  line  (5)  according  to  conditions 
(Bl),  (B2)  and  (B3)  (chapter  6.5.2).  Function  end  (x)  in 
procedure  LISU  (figure  6.5)  is  defined  by 

end(x) :=  (h ' (x)  =0)  , 

and  heuristic  h'  of  a  pixel  x  =  (col,  row)  is  computed  as 
follows 
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h' (x) :=  min 


The  cost  function  c  has  the  form  as  described  in  chapter 
6.5.2. 


(col-gc) 2+ (row-gr) 2  I  (gc,gr)  €  goal 


} 


6.5.4  Performance  of  sequential  line  detection 


Sequential  line  detection  has  its  origin  in  a  general 
problem  solving  algorithm,  which  searches  in  a  predefined 
search  space  of  possible  solutions  for  the  best  under  a  cer¬ 
tain  valuation  criterion.  How  fast  algorithm  finds  the 
result  depends  essentially  on  the  choice  of  heurictic  h1 
or  as  in  our  case  -  of  the  deduced  cost  function.  If  the  per¬ 
fect  heuristic  h  were  known,  that  is  if  condition  (B)  in  chap¬ 
ter  6.5.2  takes  equality  on  the  right  side,  then  the  algorithm 
would  find  the  optimal  solution  (xo,xlf . . . ,xn)  exactly  in  n  steps. 

A  disadvantage  of  these  algorithms  is  the  enormous 
computational  effort  if  the  line  does  not  lie  in  the  expected 
region.  Then  almost  the  whole  search  space  must  be  inspected. 

It  is  therefore  of  great  interest  to  keep  the  search  space  X 
as  small  as  possible. 


6.6  VERIFICATION  OF  FEATURES 


Features  or  objects  appearing  in  maps  and  images  are 
determined  by  several  characteristics  (form,  uniform  grey 
value) .  In  step  1  of  the  recognition  process  some  of  them  were  ' 
used  to  extract  a  feature  (binary  matrix)  from  the  image. 
Verification  is  then  the  second  step  which  tests  if  other  charac¬ 
teristics  describing  the  object  have  been  violated  in  step  1 
or  not. 

Since  three  of  our  recognition  procedures  (THRESH, ADAPT, 
LINDET)  are  adapting  the  form  of  the  object  to  yield  the  best 
attainable  uniform  grey  scale  distribution,  the  resulting 
shape  of  the  feature  must  be  compared  to  the  original  one 
from  the  map.  Figure  6.7  shows  in  lines  (7)  until  (10)  the 
general  method  of  operation.  Measures  of  the  object  and  the 
pattern  are  computed  in  the  appropriate  coordinate  system 
in  order  to  construct  a  linear  transformation  function  T 
which  superimposes  the  transformed  object  over  the  pattern. 

Then  the  shape  of  the  two  features  can  be  compared  and  measured 
by  a  correlation  function. 

Some  of  these  measures,  their  computation  and  the  recon¬ 
struction  of  a  linear  transformation  will  be  described  in 
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chapter  6.6.1.  Correlation  functions  to  measure  similarity 
of  the  two  features  are  presented  in  chapter  6.6.2. 


(1)  procedure  VERIFICATION  (object,  image,  pattern) ; 

(2)  begin 

(3)  if  search  method  (object) =  1  then 

(4)  VERIFICATION :=  coco (image,  pattern) 

(5)  else 

(6)  begin 

(7)  objectmeasures :=  measures  (object); 

(8)  featuremeasures  :  =  measures  (pattern) ; 

(9)  BESTOVERLAY  (objectmeasures,  featuremeasures,  T) ; 

(10)  VERIFICATION : =  correlation  (T (object) , pattern) ; 

(11)  end; 

(12)  end« 


Figure  6.7:  Computation  of  a  measure  for  verification. 


The  first  recognition  procedure  SHIFT  does  not  vary 
the  shape  of  the  feature.  It  is  therefore  meaningless  to 
compare  the  two  (in  that  case  identical)  shapes.  We  use 
either  the  correlation  between  grey  values  of  the  picture 
and  the  binary  pattern  that  has  already  been  calculated  in 
procedure  SHIFT  or  we  correlate  the  pattern  with  another 
channel  of  the  multiple  (multispectral)  image  or  a  pre- 
processed  version  of  it. 


6.6.1  Some  measures  describing  patterns,  their  computation 
and  application. 

(a)  Translation  and  Center  of  Gravity 


The  measures  or  quantities  discussed  here  have  the 
purpose  to  establish  a  linear  transformation 

1  is)  ~  A  X  +  h 

from  the  map  coordinate  system  to  the  image  system.  These 
measures  are  a  derived  quantity,  obtained  from  possibly 
erroneous  pixels  that  are  defined  to  belong  to  a  feature. 
Errors  in  the  pixels  should  have  a  minimal  effect  on  the 
derived  measure . 

The  simplest  measure  we  use  is  the  center  of  gravity 
of  a  pattern.  Its  computation  for  a  closed  polygon  P 
employes  the  formula 

Sx  ( P )  :  —  i[|Z(xi+xi+1)(xi+1yi  -  xiYi+1)]  /  2Z(*i+1  y±  -xiyi+l) 
Sy  (P)  :=  ^(Yi+Yi  +  iMKi+iyi  “  xiyi+l)]  7  ^  (xi+lyi  "  xiYi+l) 
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where  (x±,  y.),  i  =  1,  .  ..n  are  the  points  of  the  polygon  P, 
and  (xn+^,  y^+i)  is  identical  to  (x^  yx)  . 

If  a  region  R  is  described  by  two  polygons  P^  and  P2, 

Pj^  enclosed  in  P2  (see  figure  6.8),  then  the 


Figure  6.8:  The  shape  of  region  R  containing  an  island 
is  described  by  the  polygon  Pj^  and  P2> 


center  of  gravity  s (R)  of  R  is  computed  following 


S(P,)*F(P.)  -  S(P.)*F(P.) 

S(R)  :=  - * - * - i - 

F(P2)  -  F  (P1) 

with  F(P. ) ,  F(P2)  denote  the  areas  of  polygons  P^  and  P2, 
respectively . 

The  center  of  gravity  of  a  binary  matrix  M  computes 
as  follows 

S(M)  :=  (  S  +  M(*))  /  F(M) 

xeWindow{M) 


where  F (M)  is  here  the  number  of  pixels  in  the  pattern 
(M(x)  =  l)  . 

With  the  two  centers  of  gravity,  namely  of  region  R  in 
the  map  £}.(£) ,  and  of  pattern  M  in  the  image,  S. CM)  ,  the  shift 
vector  of  the  original  transformation  TQ 

IqCx)  -  Aq  •  x  +  ^ 


can  be  corrected  to  be 

b.T  :=  S  (M)  -  A  •  S  (R)  . 

The  new  transformation  T„  can  be  formulated  as  follows 

N 

Tn(£):=  b0  •  &  -  £(*>  )  +  i(M)  • 
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(b)  Rotation  and  scale  with  l.ines. 

Other  measures  can  be  calculated  from  lines,  two  refe¬ 
rence  points  £,  £  in  image  coordinates  and  £'  and  £'  in 
map  coordinates  are  used  to  correct  the  original  transforma¬ 
tion 

v*>  -  A,  .  a  +  So 

to 

*  *  +  *N 

with 


The  correction  matrice  D  can  be  calculated  from 

2-  £=^02)  -  In(£,)=  -  5o(£')  ). 

Having  the  new  £N,  JaN  computes  as  for  translation: 

=  £  -  An  ■  £' 


6.6.2  Measures  of  correlation  between  object  and  pattern 

The  measures  mentioned  here  are  used  to  determine  how 
well  the  shape  of  the  (best-) transformed  object  corresponds 
to  the  shape  of  the  pattern  in  the  image.  We  assume  that 
the  object  is  transformed  linearly  so  that  the  registration 
of  object  and  pattern  is  optimal. 

Then  the  polygonal  representation  is  transformed  to  the 
raster  of  the  pattern.  A  measure  can  be  obtained  using  the 
binary  coefficient  of  correlation  from  equ.(6.9),  "bicoco", 
or  the  binary  cross-correlation  from  equ.  (6.10),  "bicross", 
of  the  two  binary  matrices  (see  chapter  6.1). 

Another  method  we  used  earlier  was  based  on  the  number 
"differ"  of  pixels  which  do  not  appear  in  both  matrices  M.. 
and  M2 : 

F  (  J  M.  -  M,  |  ) 

differ:*  - — - * - 

F  (  Mx  ) 

Tests  showed  that  this  approach  is  very  sensitive  to  the 
total  surface  of  Mo.  The  best  results  were  obtained  from  the 
cross-correlation  "bicross"  which  is  simple  to  calculate 
and  less  dependent  from  the  total  area  of  a  pattern  or  object. 
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7.  CONTROL  POINTS  AND  RECTIFICATION 


The  original  problem  of  registering  a  map  with  the  cor¬ 
responding  satellite  image  was  partitioned  in  our  approach 
to  search  first  for  registration  of  separate  identifiable 
objects.  From  each  such  object  we  extract  one  or  more  pairs 
of  corresponding  points  in  map  and  image,  which  we  call 
control  points.  During  the  process  of  recognizing  objects  in 
the  image  these  control  points  are  collected.  If  the  set  of 
control  points  covers  the  overlay  region  with  sufficient  den¬ 
sity  and  an  even  distribution,  it  is  used  to  geometrically 
transform  the  image  to  the  map  coordinate  system. 

Procedure  CONTROLP  (fig.  7.1)  summarizes  the  selection 
of  control  points  from  an  identified  object  and  their  collec¬ 
tion  in  a  file.  Regions  and  lines  are  handled  separately  be¬ 
cause  of  the  difference  in  structure.  From  a  region  we  ex¬ 
tract  one  controlpoint,  the  center  of  gravity.  The  degenerated 
case  of  a  circular  region  demonstrates  why  we  restrict  ourselves 
to  only  one  point.  From  lines  we  extract  at  least  two  points. 


(1)  procedure  CCNTRDLP  (object,  pattern,  oontrolpoints) 

(2)  begin 

(3)  if  search  method  (object)  =  4  then 

(4)  begin 

0  Lines  0; 

( 5 )  oontrolpoints :  ==oontrolpoints+  &.  (object)  ,£,  (pattern) ) 

+(£,  (object)  (pattern) )  ,* 

(6)  end  else  2 

(7)  begin 

0  regions  0; 

(8)  oontrolpoints :=ocntrolpoints+(§ (object)  ,5. (pattern) ) ; 

(9)  end 

(10)  end. 


Figure  7.1:  Definition  of  control  points 

We  could  use  the  endpoints,  but  they  are  affected  with  errors 
because  lines  are  often  segments  of  longer  lines  for  which 
the  starting  point  in  the  image  is  not  exactly  defined.  To 
eliminate  this  disadvantage  and  source  of  errors  the  line  is 
subdivided  into  three  segments  of  equal  length.  The  intersec¬ 
tion  points  I-,  I2  of  adjacent  segments  are  then  used  as  con¬ 
trol  points. 

The  registration  process  REGISTER  is  followed  by  a  geometric 
transformation  RECTIFY  (figure  7.2)  of  the  image  into  the  coor¬ 
dinate  system  of  the  map.  The  transformation  is  done  in  a  stand¬ 
ard  procedure  with  the  aid  of  the  control  points.  These  are  used 
to  define  a  regular  grid  of  anchor  points  in  the  rectified 
and  in  the  distorted  image.  The  grid  is  stored  as  a  matrix  of 
point  coordinates  in  the  original  image.  Rectification  is  done  by 
bi-linear  resampling  within  each  grid  mesh  as  defined  by  the  anchor 
points.  While  the  geometric  location  of  the  center  of  a  new 
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(1)  procedure  RECTIFY  (  old  image,  controlpoints ,  new  image) 

(2)  begin 

(3)  determine  rectangular  transformation  grid  (controlpoints) ; 

(4)  for  all  rectangles  of  grid  do 

(5)  begin 

(6)  for  all  pixel  of  rectangle  do 

(7)  begin 

(8)  interpolate  position  (pixel)  linearly  in  old  image; 

(9)  new  image (pixel) :=  cldimege  (position (pixel) ) ? 

(10)  end 

(11)  end 

(12)  end 


Figure  7.2:  Geometrical  rectification  of  an  image 


pixel  in  the  distorted  image  is  computed  by  bi-linear  trans¬ 
formation,  its  grey  value  is  assigned  as  that  of  the  next  neigh¬ 
bour.  The  entire  rectification  is  a  well-documented  standard 
routine  (e.g.  KONECNY,  19.75)  and  was  programmed  in  our  case 
by  WIESEL  (1977) . 
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8.  RESULTS  AND  EXPERIENCES 

8.1  INTRODUCTION 


Numerical  results  and  experiences  were  obtained  using 
3  different  images: 

(a)  A  synthetic  image  (denoted  MERK) . 

(b)  A  section  of  a  multispectral  Landsat  image  of  an  area 
in  Southern  Germany  with  lakes  (denoted  WALCH) . 

(c)  A  section  of  a  Landsat  image  in  Southern  Austria  of  an 
urban  area  (denoted  GRAZ) . 

Tests  and  results  will  be  described  in  this  chapter? 
they  will  demonstrate  in  some  detail  the  basic  concept  of 
ARSIM.  In  each  of  the  three  test  images  several  specific 
aspects  of  the  ARSIM  concept  have  been  evaluated. 

Experiments  with  the  synthetic  image  MERK  serve  to 
prove  the  basic  validity  of  the  procedures,  using  a  well 
controlled  ideal  data  set. 

Test  area  WALCH  was  used  originally  for  understanding 
the  problem  of  feature  recognition  in  satellite  images.  It 
served  to  develop  the  rudimentary  algorithms  and  to  select 
methods  which  promised  some  success. 

Test  area  GRAZ  finally  had  the  function  of  verifying 
the  developed  methods  with  real  data.  This  is  the  material 
to  permit  identification  of  current  limits  and  needed 
improvements  in  the  experimental  system.  The  fundamental 
experiences  can  be  summarized  as  follows: 

-  ARSIM  is  a  valid  concept  to  automatically  merge  a  map 
with  an  image. 

The  task  of  automatically  merging  maps  and  satellite 
images  is  a  formidable  one.  There  is  no  limit  to  the 
number  of  possible  pattern  recognition  techniques  that 
can  be  used. 

Successful  feature  recognition  depends  on  specific 
preprocessing  to  prepare  the  feature  for  a  recognition 
algorithm.  Preprocessing  must  be  turned  to  both  the 
image  and  to  the  object  to  be  recognized. 

Successful  feature  recognition  also  depends  on  a  map 
data  base  of  an  appropriate  organisation  and  content 
such  that  a  relation  to  images  can  be  established. 

The  success  of  the  approach  depends  on  the  implementation 
of  experiences  into  the  selection  of  preprocessing  and 
recognition  algorithms.  A  "learning  system"  is  required. 

ARSIM  must  be  organized  in  a  highly  modular  form  so  that 
many  different  algorithms  can  be  implemented  as  new 
experiences  are  gained. 
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8.2  EXPERIENCES  WITH  THE  SYNTHETIC  IMAGE 


The  purpose  of  working  with  a  synthetic  image  and 
related  data  base  is  to  have  control  over  perturbation 
factors.  Therefore  the  results  can  be  more  easily  analyzed 
and  the  limitations  of  the  concept  can  be  understood. 

8.2.1  The  synthetic  image  MERK 

Figure  8.1  shows  in  the  form  of  a  printer  output  the 
configuration  of  the  five  objects  RECHT,  DREI,  KARO,  PARA, 

KREIS  in  the  image.  Different  geometric  forms  of  objects 
have  been  used:  RECHT  is  a  rectangle,  DREI  a  triangle,  KARO 
a  quadrangle  with  a  rectangular  hole,  PARA  a  long  parallelogram 
and  KREIS  is  a  circle.  Pixels  of  the  objects  were  set  to  grey 
value  0  (black)  and  the  background  to  255  (white) . 

A  map  data  base  was  created  showing  the  same  features 
with  identical  geometry.  Figures  8.2  and  8.3  illustrate  the 
map. 


Project  funds  did  not  permit  us  to  fully  exploit  the 
possibilities  of  the  synthetic  images.  Instead  of  a  wide 
variety  of  perturbed  versions,  we  could  only  test  the  influence 
of  grey  level  noise  generating  a  noisy  picture  as  shown  in 
figure  8.4. 


8.2.2  A  projection  of  the  map  onto  the  image 

To  simulate  the  geometric  conditions  from  an  approximately 
rectified  satellite  image  the  map  data  were  geometrically 
deformed  using  the  following  linear  transformation  (referenced 
as  projection  (8.1)  in  the  following  chapters) : 


Tx(x,y)  =  -0.035  *  x  +  1.033  *y  +  3.500 
Ty(x,y)  =  1.058  *x  +  0.040  *y  -  5.000 


(8.1) 


Figure  8.5  shows  the  projection  of  the  five  objects  into  the 
image.  Table  8.  1  enumerates  the  objects  with  their  surface 
in  the  map  and  after  transformation  (8.1)  including  a 
discretization  error  and  the  deviation  of  the  center  of 
gravity  from  the  ideal.  These  deviations  are  measured  with 
the  Euclidean  distance.  The  assumption  is  met  of  a  maximum 
shift  of  10  pixels  in  each  coordinate  direction,  except  for 
object  DREI  which  is  used  to  study  a  more  extreme  case  of 
geometric  deformation. 


Feature 

Surface 

orig. 

Surface 

project. 

RECHT 

36.1 

41.9 

5.099 

DREI 

20.0 

24.6 

11.700 

KARO 

48-4 

44.6 

7.211 

PARA 

74.1 

87.5 

4.123 

KREIS 

48.5 

55.8 

12.042 

Table  8.1:  Characteristics  of  the  projection  in  pixel 
units . 


. 4 1  Synthetic  image  MERK  with  noise. 


8.2.3  Results  with  procedure  SHIFT 


Procedure  SHIFT  was  applied  as  discribed  in  chapter  6.2. 
For  all  objects  a  position  was  found  of  maximum  correlation 
between  the  binary  map-derived  image  and  the  window  in  MERK. 

A  section  of  the  matrix  containing  the  positional 
correlation  measures  in  percent  is  shown  in  figure  8.6.  The 
correlation  value  grows  from  14.2  X  at  the  original  position 
until  88.6  X  at  the  maximum  (the  correlation  measure  was 
explained  in  equ.  6.7).  100  X  is  not  reached  because  procedure 
SHIFT  takes  into  consideration  only  the  shift  parameters  of 
the  projection,  however,  not  rotation  and  scale  factors  which 
cause  the  deviation  of  the  correlation  value  from  the  ideal 
value  of  100  X. 

Procedure  SHIFT  was  applied  also  to  the  noisy  version  of 
image  MERK.  The  maximum  correlation  was  practically  the  same 
or  even  higher  in  some  instances  in  the  noisy  image.  This  can 
be  explained  by  the  smoothing  effect  of  the  noise. 


79.6  ^  85.5 


85.8  I  88.6 


77.21  82.3 


76.5 


82.3 


74.1 


65.5 


67.9 


65.1 


57.5 


52.2 


54.2 


51.9 


49.3 


42.4 


41.2 


39.3 


37.1 


34.9 


28.7 


27.4 


25.6 


14.8 


IE2P 

14.1 

5.9 

■SBfKZH 

■B 

19.2 

15.4 

9.5 

Figure  8.6:  Correlation  measures  for  object  KARO  in  %,  and 
path  of  increasing  correlation  from  a  starting  value  to  the 
maximum. 

This  noise  can  compensate  for  the  effect  of  the  rotation  error. 

Figure  8.7  compares  the  results  from  both  images  with  the 
aid  of  error  vectors  of  an  overdetermined  transformation.  The 
control  points  used  for  the  transformation  are  the  centers  of 
gravity  of  the  features  in  the  map  and  the  shifted  features  in 
the  projected  version  of  the  map. 

To  qualify  the  method  a  standard  deviation  for  coordinate 
errors  and  point  errors  (vector  lengths)  are  calculated  using: 

— 1 


-  JzZ 

»  i=l 


/ (n-3) 


where  d,  are  the  individual  residuals  after  the  transformation, 
n  is  the  number  of  points;  (n-3)  is  the  degree  of  freedom 
(overdetermination)  of  the  transformation. 
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Clean 

Feature 

mama 

HEQMKnMi 

Length 

RECHT 

DREI 

KARO 

PARA 

KREIS 

(  61.500  ,  61.520) 

(  58.667  *  183.667) 
(124.500  , 124.949) 
(186.500  ,  79.500) 
(187.000  , 187.000) 

(  0.692  ,  -0.157) 
(-0.760  ,  0.347) 

(0.178  ,  -0.400) 
(-0.887  ,  0.404) 

(  0.777  ,  -0.194) 

0.709 

0.835 

0.437 

0.974 

0.800 

87.230 

72.783 

88.575 

88.969 

88.765 

Stand,  deviation  (1.113  ,  0.503) 

1.221 

85.264 

Noise 

RECHT 

DREI 

KARO 

PARA 

KREIS 

(  61.500  ,  61.500) 

(  58.667  ,183.667) 
(124.500  ,124.949) 
(186.500  ,  79.500) 
(187.000  , 187.000) 

(  0.692  ,  -0.073) 
(-0.760  ,  0.220) 

(  0.178  ,  -0.307) 
(-0.877  ,  0.255) 

(0.777  ,  -0.095) 

0.696 

0.791 

0.355 

0.923 

0.782 

87.160 

72.855 

88.592 

88.878 

88.777 

Stand,  deviation  (  1.113  ,  0.333) 

1.162 

85.252 

Figure  8.7:  Table  of  results  from  procedure  SHIFT. 


The  comparison  of  the  clean-  and  noisy  data  shows  that 
maximum  correlation  appears  always  at  the  same  position 
except  for  object  DREI.  Object  DREI  finds  maximum  correlation 
at  the  border  of  the  area  of  interest,  consisting  of  21  x  21 
pixels  in  a  square  around  the  original  position.  The  absolute 
maximum  would  lie  outside  this  window.  In  the  clean  version 
two  adjacent  positions  get  the  same  maximum  value  while  in 
the  other  version  noise  influences  the  correlation.  Hence 
the  maximum  is  found  in  only  one  position.  This  is 
reflected  in  the  change  of  one  coordinate  of  the  control  points 
of  DREI  and  explains  why  the  error  vectors  are  not  identical 
in  the  two  versions. 


8.2.4  Results  with  procedure  THRESH 

In  the  same  way  as  procedure  SHIFT,  procedure  THRESH  was 
applied  to  the  two  versions  of  image  MERK  with  projection  (8.1) 

The  parameters  chosen  for  this  example  were  distance  D  * 
10  and  factor  p  =  1.0  (compare  chapter  6.3).  Two  restrictions 
to  the  threshold  bounds  IPU  and  IPO  were  necessary: 

(a)  The  lower  bound  IPU  was  fixed  to  0  and 

(b)  "bad"  grey  values  were  not  allowed  by  setting  IPO  <  255. 

Tests  have  shown  that  no  reasonable  results  can  be 
obtained  without  these  two  restrictions  since  the  overlap 
region  between  projection  and  image  feature  is  so  small.  An 
alternative  has  been  to  choose  factor  p  so  small  that  p*F (M) 
is  less  than  the  surface  of  the  overlap.  This  works  for  the 
clean  version  but  not  with  noise. 
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The  two  conditions  above  seem  to  be  reasonable  because  they 
can  be  translated  into  conditions  for  preprocessing  which 
underlying  satellite  images  should  satisfy: 

(a)  The  grey  values  of  the  preprocessed  image  should  be  low 
for  pixels  which  belong  to  the  feature  and  high  for  the 
rest. 

(b)  Pixels  with  grey  value  255  are  excluded  a  priori. 

Figure  8.8  summarizes  the  results  of  THRESH. 


C 

Lea 

n 

Feature 

mamm 

— 

■nM 

Length 

Correl . 
[X] 

RECHT 

(  61.500  ,  61.500) 

(  0.223  , 

-0.224) 

0.316 

98.533 

DREI 

(  58.667  , 183.667) 

(-0.361  , 

0.248) 

0.438 

85.217 

KARO 

(124.500  ,124.949) 

(  0.298  , 

-0.062) 

0.304 

92.823 

PARA 

(186.500  ,  79.500) 

(-0.421  , 

0.289) 

0.510 

93.439 

KREIS 

(187.000  ,187.000 

(  0.261  , 

-0.251) 

0.363 

94.913 

Stand,  deviation 

(  0.507  , 

0.362) 

0.623 

92.805 

N  < 

3  i  s 

e 

RECHT 

(  61.500  ,  61.500)' 

(  0.662  , 

-0.178) 

0.685 

89.472 

DREI 

(  58.667  ,183.667) 

(-0.721  , 

0.268) 

0.770 

43.128 

KARO 

(124.500  ,124.949) 

(  0.159  , 

-0.195) 

0.252 

65.279 

PARA 

(186.500  ,  79.500) 

(-0.842  , 

0.312) 

0.898 

86.673 

KREIS 

(187.000  ,187.000) 

(  0.743  , 

-0.207) 

0.771 

77.387 

Stand.,  deviation 

(  1.059  , 

0.375) 

1.124 

72.388 

Figure  8.8:  Table  of  results  from  procedure  THRESH. 


Binary  correlation  is  used  here  for  verification.  100  *  is 
not  reached  for  the  same  reason  as  for  SHIFT.  However,  better 
values  for  the  clean  image  show  that  thresholding  takes  care 
of  feature  rotation  and  scale  in  the  image,  but  verification 
considers  only  the  shift  of  the  projection. 

Results  from  the  noisy  image  demonstrate  the  strong  noise 
dependency  of  thresholding,  on  the  average  correlation  is 
20  X  smaller  than  before.  Moreover  object  DREI  could  not  be 
accepted  as  recognized  due  to  a  correlation  value  of  43x. 

Figure  8.9  explains  why  DREI  has  not  been  recognized  automatically 
There  is  a  very  small  overlap  between  the  image  feature  and  the 
projected  map  element  for  the  determination  of  the_ bounds. 

Hence  very  much  noise  is  included.  The  center  of  gravity  of 
all  threshold  pixels  is  therefore  displaced  so  that  verification 
cannot  state  a  similarity. 

8.2.5  Results  of  procedures  LISU  and  LINDET 

The  left  hand  side  of  the  parallelogram  PARA  was  used 
to  test  line  following.  The  edge  operator  of  Roberts  (chapter 
5.2.4)  was  applied  to  the  image  (noisy  and  clean  version:  see 
f igure  8.10) . 
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Figure  8.9x  Recognition  of  DREI  by  thresholding. 

(c)  (e)  binary  masked  printout 

(b)  (d)  (f)  complemented  binary 

original  threshold  verified 

projection (8. 1)  pattern  projection 


These  edge  images  were  used  for  line  detection.  The  starting 
mask  was  created  and  the  histogram  provided  the  following 
parameters  of  the  cost-function  (see  chapter  6.5.3): 


Clean  version 

Noisy  version 

a  =  256.0 

b  =  1.0 

n  =  0.78 

a  =  188.52 

b  =  -  0.7364 

n  =  0.85 

Table  8.2:  Parameters  of  the  costf unction  for  line  following 


Starting  and  end  point  S  and  E  are  projected  into  the  image 
(S  , E  )  and  improved  (S_,E  )  by  the  procedure  as  illustrated 
inpfi§ure  8.11:  1  A 
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Figure  8.11:  Improvement  of  starting  and  end  point  of  the 
line. 

(a)  (cost  Sj)  <  1, (cost  Ej)  <  1 

(b)  Sj  and  Ej  are  contained  in  the  starting  mask 

(c)  Among  pixels  satisfying  conditions  (a)  and  (b)  S_  and  E 

have  the  maximum  distance  to  and  E^. 
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Black  pixels  in  figure  8.11  are  pixels  satisfying  the  first 
two  conditions.  The  arrow  points  to  the  improved  extreme 
points  corresponding  to  condition  3. 

Figure  8.12  shows  the  results  of  procedures  LISU  and 
LINDET  within  the  corresponding  search  area. 

For  verification  purposes  the  corresponding  lines  in 
the  map  and  image  were  subdivided  into  6  parts  of  equal  length. 
Subdivision  points  2  and  4  are  used  to  calculate  a  new 
transformation  of  best  overlay  according  to 
chapter  6.6.1(b).  All  5  subdivision  points  of  the  map  are 
transformed  and  difference  vectors  calculated  to  the  5 
corresponding  points  of  the  image  line.  The  results  are  shown 
in  figure  8.13.  No  errors  appear  in  the  clean  version  because 
the  line  is  straight.  In  the  noisy  picture  some  pixels  do  not 
lie  on  a  straight  line.  Starting  and  end  points  are  not 
included  in  the  verification  process  because  they  are  not  a 
result  of  line  following. 


Point  differences  of  the  linear  transformation  adapted  to  points  2  and  4 

Subdivision 

points 

Error  (clean 
Image) 
Ax  Av 

length 

Error  (noisy  image) 

Ax  Ay 

length 

1 

(.0  ,  .0) 

.0 

(  .207  ,  -.207) 

.293 

2 

(.0  ,  .0) 

.0 

(  .0  ,  .0  ) 

.0 

3 

(-0  ,  .0) 

.0 

(-.146  ,  .232) 

.275 

4 

(.0  ,  .0) 

.0 

(  .0  ,  .0  ) 

.0 

5 

(.0  ,  .0) 

.0 

(-.207  ,  .207) 

.293 

Stand. dev. 

(.0  ,  .0) 

.0 

(  .231  ,  .264) 

.352 

Figure  8.13:  Verification  of  the  line  in  MERK. 


8.2.6  Results  with  procedure  ADAPT 


Considerations  of  computing  time  caused  us  to  use  a  new 
approximate  projection  which  is  improved  and  brings  the 
objects  nearer  to  the  image  features.  This  projection  (8.2) 
has  the  following  parameters 


T  (x,y)  =  -0.01  *  x  +  1.02  *  y  -  0.50 

A 

Ty(x,y)  =  1.02  *  x  -  1.00 


(8.2) 


and  results  in  areas  for  each  feature  as  shown  in  table  8.3. 
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Feature 

Area 

orig. 

Area 

proj. 

Differences  in 
position (pixels ) 

RECHT 

36.1 

40.0 

0.000 

DREI 

20.0 

23.2 

2.968 

KARO 

48.4 

41.0 

1.759 

PARA 

74.1 

83.2 

2.570 

KREIS 

48.5  . 

53.9 

2.998 

Table  8.3:  Characteristics  of  projection  2  dn  pixel  units 


Procedure  ADAPT  is  demonstrated  by  means  of  object  KARO 
(figure  8.14).  Four  iterations  have  been  computed.  The 
changes  of  the  shapes  are  illustrated  in  figure  8.15.  This 
figure  presents  the  surface  of  the  feature  "binadapt"  before 
each  iteration  step;  the  number  of  pixels  that  are  added  (+) 
and  eliminated  (-)  at  each  step;  and  finally  the  grey  value 
bounds  and  the  number  of  pixels  contained  in  this  interval 
in  binadapt. 


Iteration 

Starting 

Area 

+ 

- 

Grey  value 
limits :  area 

1 

410 

32 

48 

(0,254) :334 

2 

394 

20 

17 

(0,254) : 366 

3 

397 

0 

11 

(0,254) : 386 

4 

386 

0 

00 

(0,0  ) : 386 

386 

Figure  8.15:  Shape  changes  while  iterating  ADAPT  on 
oject  KARO. 


During  the  last  iteration  step  the  correct  interval  (0,0)  is 
found.  Verification,  however,  is  less  than  100  X  because  of 
orientation  and  scale  differences  in  projection  (8.2)  (see 
figure  8. 14e, f ) . 

A  summary  of  results  from  procedure  ADAPT  is  presented 
in  figure  8.16. 
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C 

lea 

n 

Feature 

Position  in  the  image 

Error 

vector 

Length 

Correl. 

[*] 

x  y 

AX 

Ay 

RECHT 

(  61.500  ,  61.500) 

(  0.134  , 

-0.010) 

0.134 

100.000 

DREI 

(  58.667  , 183.667) 

(-0.133  , 

0.015) 

0.134 

99.349 

KARO 

(124.500  , 124.949) 

(  0.006  « 

-0.011) 

0.012 

95.867 

PARA 

(186.500  .  79.500) 

(-0.155  , 

0.017) 

0.156 

98.533 

KREIS 

(187.000  ,  187.000) 

(  0.149  , 

-0.011) 

0.149 

98.094 

Stand.,  deviation 

(  0.202  , 

0.021) 

0.203 

98.369 

N 

o  i  s 

e 

RECHT 

(  61.500  ,  61.500) 

(  0.171  , 

0.044) 

0.176 

100 . 000 

DREI 

(  58.667  , 183.667) 

(-0.127  , 

-0.058) 

0.139 

83.604 

KARO 

(124.500  . 124.949) 

(-0.082  , 

0.032) 

0.088 

83.733 

PARA 

(186.500  ,  79.500) 

(-0.148  , 

-0.068) 

0.163 

86.561 

KREIS 

(187.000  ,  187.000) 

(  0.186  , 

0.050) 

0.193 

87.343 

Stand,  deviation 

(  0.233  , 

0.082) 

0.247 

88.248 

Figure  8.16:  Results  from  procedure  ADAPT. 


For  object  RECHT  projection  (8.2)  results  in  an  errorfree 
overlap  which  is  detected  by  ADAPT  in  the  first  iteration. 
Noise  decreases  the  correlation  by  about  10  £.  The  effect  of 
noise  in  the  control  points  is  very  small. 


8.2.7  Comparison  of  recognition  procedures 

After  the  determination  of  control  points  for  every 
object  with  every  of  the  investigated  methods  an  overdetermined 
6-parameter  linear  transformation  is  computed  with  residuals 
to  compare  the  results.  Figure  8.17  presents  the  results  from 
the  clean  and  noisy  image. 


c  1  e 

a  n 

n  o  i  s 

e 

method 

feature 

error  vector 

length 

error  vector 

length 

SHIFT 

RECHT 

(  0.450  .  -0.657) 

0.796 

(  0.905  ,  -0.567) 

1.068 

DREI 

( -1.712  ,  1.235) 

2.111 

(  0.117  ,  1.735) 

1.739 

KARO 

(  0.083  ,  -0.517) 

0.524 

(  0.808  ,  0.175) 

0.827 

PARA 

(  -0.229  ,  -0.550) 

0.595 

(  -0.434  ,  -0.074) 

0.440 

KREIS 

(  0.833  ,  0.058) 

0.835 

(  1.816  ,  1.340) 

2.257 

THRESH 

RECHT 

(  0.470  ,  -0.181) 

0. 504 

(  0.089  ,  -0.275) 

0.289 

DREI 

(  -0.655  ,  0.344) 

0.740 

(  -3.396  ,  -0.870) 

3.506 

KARO 

(  0.323  .  -0.294) 

0.437 

(  -0.982  ,  -1.008) 

1.407 

PARA 

(  -0.139  ,  -0.250) 

0.286 

(  -0.706  ,  -0.278) 

0.759 

KREIS 

(  0.070  ,  -0.757) 

0.760 

(  -0.947  ,  -1.723) 

1.966 

Figure  8.17:  Overall  results  with  synthetic  image  MERK  (cont'd). 
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ADAPT 

RECHT 

DREI 

KARO 

PARA 

KREIS 

(  0.470  ,  -0.181) 

(  0.665  ,  -0.138) 

(  0.383  ,  -0.223) 
(-0.139  .  -0.250) 

(  0.562  ,  -0.265) 

0.504 

0.679 

0.443 

0.286 

0.621 

(  0.925  ,  -0.091) 

(  1.131  ,  0.535) 

(  0.780  ,  -0.284) 

(  0.382  ,  -1.184) 

(  1.152  /  -0.452) 

0.930 

1.251 

0.830 

1.244 

1.237 

LINDET 

PUNKT 2 
PUNKT 4 

(-0.383  ,  0.789) 

(-1.051  ,  1.835) 

0.877 

2.115 

C-0.674  t  1.019) 
(-0.966  ,  2.001) 

1.222 

2.222 

Stand. deviation 

(  0.710,  0.738) 

1.024 

(  1.322  /  1.105) 

1.723 

Figure  8.17:  Comparison  of  recognition  methods  by  error  vectors. 


For  lines,  the  subdivision  points  2  and  4  are  used  as 
control  points.  In  the  clean  version  there  are  two  points 
with  error  vector  lengths  greater  than  1:  DREI  with  SHIFT 
and  point  4  with  LINDET.  DREI  was  the  object  from  which  the 
projection  (8.1)  was  more  than  10  pixels  apart  from  the  image 
feature.  The  error  of  the  line  subdivision  point  4  has  its 
explanation  in  the  fact  that  the  end  point  E_  in  the  image 
lies  on  the  short  side  of  the  parallelogram  XPARA  while  the 
starting  point  was  found  correctly.  So  the  line  in  the  image 
must  be  longer  than  the  real  projection  of  the  map.  The  effect 
of  that  is  that  the  control  points  are  shifted  towards  the 
end  point. 

From  the  point  of  view  of  noise  all  errors  of  all 
procedures  become  greater  but  THRESH  seems  to  be  the  worst. 
Errors  of  LINDET  are  of  the  same  type  as  in  the  clean  image. 
Errors  of  ADAPT  seem  to  be  the  most  homogeneous  (0.2  to  0.7 
with  clean  data,  0.8  to  1.2  with  noise)  with  little  variations 
while  the  other  methods  vary  more  in  their  performance. 


Figure  8.18  summarizes  the  standard  deviation  of  all 
methods  and  the  mean  of  the  correlation  measures  used  for 
verification. 


standard  deviations 

n 

method 

noise 

0x 

mm 

°P 

mean 

correlatior 

SHIFT 

no 

1.113 

0.503 

1.221 

85.264 

X 

yes 

1.113 

0.333 

1.162 

85.252 

X 

THRESH 

no 

0.507 

0.362 

0.623 

92.805 

X 

yes 

1.059 

0.375 

1.124 

72.388 

X 

ADAPT 

no 

0.202 

0.021 

0.203 

98.369 

X 

yes 

0.233 

0.082 

0.247 

88.248 

X 

LINDET 

no 

0.000 

0.000 

0.000 

100.000 

X 

yes 

0.231 

0.264 

0.352 

98.278 

X 

Figure  8.18:  Comparison  of  recognition  methods  by  statistical 
measures . 


-121- 


The  correlation  measure  licorr  of  the  line  is  calculated 
from  the  error  vector  lengths  liwiththe  following  formula: 


licorr 


100 


(,.  (  tuU«)) 


where  n  is  the  number  of  intersection  points  and  D  is  the 
maximum  shift  lengths  as  assumed. 

This  table  shows  again  that  method  THRESH  is  most 
sensitive  to  noise.  In  general  it  is  better  than  shift  and 
has  its  greatest  advantage  in  computation  speed.  Speed  was 
the  handicap  of  procedure  ADAPT  in  our  tests.  However,  the 
current  version  of  ADAPT  provides  plenty  of  possibilities 
and  parameters  to  accelerate  speed  in  an  improved  version. 


8.3  EXPERIENCES  WITH  LANDSAT  IMAGE  WALCH 


Image  WALCH  was  the  first  satellite  image  treated. 

Basic  concepts  of  the  recognition  methods  were  implemented 
and  tested  and  often  corrected.  The  purpose  was  to  obtain 
an  understanding  of  weaknesses,  not,  however,  to  compare 
the  methods. 

8.3.1  Image  WALCH  and  map 

Figure  8.19  shows  the  portion  of  a  satellite  image  from 
Southern  Germany  in  its  4  spectral  bands  MSS4,  MSS5,  MSS6 
and  MSS 7. 

This  image  was  available  in  a  NASA-format.  This  is 
significant  for  the  present  purpose,  since  NASA  performs  a 
grey  value  normalisation.  Grey  values  in  MSS4,  MSS5  and 
MSS6  range  from  0  to  127,  while  in  MSS7  it  is  between  0  and 
63. 


•  The  corresponding  map  data  bank  covers  the  area  of  the 
image.  Figure  8,20  shows  a  plot  output  of  the  map.  Indicated 
are  some  features  which  were  used  for  recognition. 

The  features  are  of  different  size,  texture  and  grey 
tone.  There  are  villages  (WA-SEE-ORT,  Schledorf ) ,  there  is 
an  island,  a  forest  (Karpfau) ,  grassland  (Jochberg,  Sachen- 
bach,  Wasuedufer,  Wasuedufer  2,  Wasuedufer  3,  Kochelnord) 
and  water  bodies  (Walchensee,  Kochelsee,  Karpfsee) . 


Lite  image 


KocRelnord 


Figure  8.20:  Plot  of  part  of  the  digital  map  data  bank  for 

the  test  area  in  Southern  Germany.  The  indicated 
features  were  used  in  the  test  of  the.  procedure . 
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Table  8.4  provides  a  list  of  objects  ordered  acc.  to 
size  .  It  could  be  the  result  of  procedure  OBJECTSELECTION 
(chapter  2.6).  The  table  also  lists  the  centers  of  gravity 
of  each  object  which  are  needed  during  the  overlap  process 
and  the  type  of  each  region  (1  =  water  body,  3  =  villages, 

4  a  grassland,  5  =  forest,  9  -  island). 


Feature 

Area  in 
pixels 

Center  of  Gravity 
x  y 

Type 

WALCHENSEE 

3589 

(1502.268,498.028) 

1 

KOCHELSEE 

1374 

(1501.242,510.522) 

1 

KOCHELNORD 

503 

(1500.632,512.897) 

4 

SCHLEDORF-UF 

339 

(1497.685,510.566) 

4 

KAKOPFWALD 

176 

(1497.627,496.215) 

5 

WA-SEE-ORT 

146 

(1496.879,498.143) 

3 

SACHENBACH 

126 

(1506.761,501.946) 

4 

JOCHBERG 

113 

(1506.301,505.659) 

4 

WALD- 10 

89 

(1496.127,494.969) 

5 

SCHLEDORF 

63 

(1497.593,512.858) 

3 

KAKOPFFREI 

32 

(1498.525,496.517) 

4 

WA-SUEDUFER 

30 

(1501.248,494.392) 

4 

WA-SUEDUFER  2 

28 

(1502.809,494.141) 

4 

KARPFSEE 

21 

(1496.203,514.633) 

1 

INSEL 

14 

(1505.347,498.391) 

9 

WA-SUEDUFER  3 

11 

(1504.169,494.278) 

4 

KARPFAU 

3 

(1496.646,514.131) 

5 

Table  8.4  :  Sorted  objectlist  for  map  WALCH. 


8.3.2  Projecting  the  map  into  the  image 

An  approximate  geometric  relationship  must  be  established 
prior  to  the  automated  recognition  process.  The  map  data  are 
therefore  projected  into  the  image  with  an  approximate  a-priori 
6-parameter  transformation.  This  approximation  may,  for  example, 
derive  from  the  satellite  orbit  data. 

The  result  of  the  transformation  is  shown  in  figure 
8.21.  The  transformation  used  the  following  parameters: 

Tx(x,y)  =  8.213448  *  x  -  2.60549  *  y  -  10685.4  {g 

Ty(x,y)  =  -1.489149  *  x  -  6.120147  *y  +  5423.8 


(a)  Walchensee 


Figure  8.21:  Projection  of  objects  WALCHENSEE,  KOCHELSEE, 

KAKOPFWALD,  SCHLEDORF-UF  and  QALD-10  into 
J -  LCH  feont'dl . 
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8.3.3  Results  from  procedure  SHIFT 

Procedure  SHIFT  was  applied  to  24  objects  of  the  map. 

13  could  be  recognized  with  a  correlation  greater  than  50  X. 
Figure  8.22  lists  the  names  of  objects,  the  channels  for 
which  maximum  correlation  was  reached,'  the  correlation 
measures  in  percent  and  error  vectors  with  their  lengths. 


Feature 

Error  vector 
x  y 

Length  of 
vector 

1 

Corre¬ 

lation 

WALCHENSEE 

(-0.690 

/ 

0.161) 

0.708 

7 

60.658 

KOCHELSEE 

(-1.107 

9 

-0.398) 

1.177 

7 

73.102 

KARPFSEE 

(  1.682 

9 

-0.162) 

1.689 

6 

55.519 

WA-SEE-ORT 

(-1.317 

9 

0.630) 

1.460 

4 

55.057 

SCHLEDORF 

(-1.328 

9 

-0.285) 

1.358 

4 

67.277 

WA-SUEDUFER 

(  0.714 

9 

0.425) 

0.831 

7 

61.290 

WA-SUEDUFER  2 

(  2.054 

9 

0.319) 

2.079 

5 

72.195 

WA-SUEDUFER  3 

(-0.943 

9 

-1.767) 

2.003 

5 

58.777 

JOCHBERG 

(-0.166 

9 

0.379) 

0.414 

4 

79.698 

KOCHELNORD 

(  0.020 

9 

0.492) 

0.492 

5 

60.084 

SACHENBACH 

(  0.572 

9 

0.490) 

0.753 

4 

69.800 

KARPFAU 

(  0.720 

9 

-0.258) 

0.765 

5 

69.479 

INSEL 

(-0.211 

9 

-0.026) 

0.213 

7 

60.999 

Stand .  deviation 

(  1.207 

9 

0.692) 

1.391 

Figure  8.22:  Summary  of  results  from  procedure  SHIFT. 


The  distribution  of  the  projected  objects  and  the 
direction  in  which  they  are  shifted  are  shown  in  figure  8.23. 

The  control  points  are  used  for  a  subsequent  image  rectification 
which  is  described  in  chapter  8.3.7. 


8.3.4  Results  from  procedure  THRESH 

Thresholding  itself  is  a  very  simple  algorithm. 
Difficulties  arise  for  the  choice  of  threshold  bounds.  How 
could  they  be  determined  automatically?  The  only  aid  Is  the 
histogram  if  no  external  information  is  available  about  grey 
values.  Even  some  rough  estimations  must  be  based  on  a 
histogram  of  grey  values  in  the  area  of  interest. 


a  no. 


Figure  8.23:  Deformation  vectors  as  obtained  from  the 
recognized  13  features  of  emthod  SHIFT. 


There  are  two  possibilities  to  vary  parameters: 

(a)  The  selection  of  a.  region  can  be  varied,  thus  the  set 
of  pixels  which  is  considered  for  the  histogram 
(histogram  mask) . 

(b)  The  algorithm  can  be  varied  to  determine  from  a  given 
histogram  the  threshold  bounds. 

The  following  histogram  masks  were  used: 

H1(D):  The  rectangular  window  enclosing  the  map  object 

as  projected  into  the  image,  plus  a  border  of  D 
pixels  width  in  each  coordinate  direction. 

H2(D):  The  projection  of  the  object  minus  a  band  of  D 

pixels  width  from  the  border. 

H3(D):  The  projection  of  the  object  plus  a  band  of  D 

pixels  width  from  the  border. 


A 


figure  8.24b. 

Figure  8.24:  Use  of  thresholding  to  collect  pixels  of  lake 

WALCHEN:  (a)  raw  thresholding  of  a  slice  of  pixels 
from  the  histogram,  (b)  noisy  pixels  are  eliminated 
using  the  map  feature  as  a  template. 
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The  sequence  thresholding-elimination  of  noisy  pixels 
via  map-guidance  was  inverted  in  the  proposed  procedure 
according  to  chapter  6.3.  The  noise  elimination  can  be 
performed  using  a  distance  transform  or  can  be  substituted 
by  the  preselection  of  the  starting  mask. 


Further  results  of  different  constellations  for  threshold 
band  determination  are  presented  in  figure  8.25  using  4  objects 


Ob  j  ect 

Thresh. 

bounds 

Correl . 

Diff.  vector 
to  shift  res. 

WALCHENSEE 

1:  HI  ( 10)  ,  Gl 

(1.26) 

54.735  X 

(-8.6  ,  -1.8) 

2:  HI ( 10) , G3 ( 1 .0) 

(0,  5) 

88.282  X 

(  0.5  ,  0.3) 

3:  H3(5)*,G3(0.86) 

(0,  5) 

86.882  X 

(  3.6  ,  0.2) 

4s  HI (9) , G3 (1. ) 

(0,  4) 

91.060  X 

(  0.7  ,  0.9) 

5:  HI (9) , G2 

(0,  8) 

7 4. 096  X 

(-1.9  ,  -1.3) 

+NE** 

89.385  X 

(-0.1  ,  -0.9) 

+NE** 

93.752  X 

(  1.0  ,  0.D 

KOCHELSEE 

Is  Hi ( 10) , G1  (5,37)  21.511  % 
2s  HI ( 10) , G3 (1 . )  (1,  5)  70.507  * 
3s  H2 (5)  *, G2  (1,  5)  67.661  % 
4s  H3 (5) *, G3 (0.77)  (1,  5)  80.866  X 


(  2.9  ,  -5.5) 
(-3.1  ,  -4.3) 
(  0.0  ,  -0.9) 
(-1.0  ,  -2.4) 


SCHLEDORF-UF*** 

Is  H3(5)*,Gl  (8,38) 

2s  H3 (5) *,G3 (0.81)  (33,45) 

KAKOPFWALD*** 

Is  H2(5)*,G1  (4,26) 

2s  H3 (5) *, Gl  (1,34) 

3s  H3(5)*,G3(0.90)  (1,  4) 


37.965  X 
51.773  X 


72.004  X 
44.982  X 
41.247  X 


Figure  8.25:  Results  from  procedure  THRESH  for  4  selected 
objects . 

*  Results  impaired  by  exceptional  effects  in  determining 
the  histogram  mask.  These  effects,  if  existing,  could 
be  eliminated  by  refined  methods. 

**  Elimination  of  noise  after  thresholding  (like  fig.  8.24). 
***  objects  not  recognized  with  SHIFT,  therefore  no  difference 
vector. 
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The  grey  tone  boundary  determination  algorithm  depends 
on  the  histogram  masks.  G1  could  be  combined  with  H2  (see 
feature  KAKOPFWALD, 1)  but  did  not  work  well  with  Hi  (WALCHEN- 
SEE, 1;  KOCHELSEE , 1 )  nor  with  H3  (KAKOPFWALD, 2)  where  the 
resulting  threshold  interval  is  too  large.  The  parameter  p 
of  G3(p)  depends  also  on  the  histogram  mask  used.  With  HI,  p 
can  be  greater  or  equal  to  1  because  the  mask 1 s  surface  is 
much  greater  than  the  surface  of  the  projection; with  H2  a 
combination  with  G3(p)  is  not  meaningful  because  the  surface 
reduction  by  H2  eliminates  already  most  of  the  noise  and  in 
this  case  p  depends  also  on  the  shape  of  the  object,  so  that  no 
proportional  relationship  can  be  established.  Good  results 
were  obtained  from  G3(p)  with  H3.  Parameter  p  was  chosen  to 
be  less  than  1.0  to  compensate  discretisation  errors  of  the 
object  projection. 

From  all  these  experiences  follows  that  the  version  of 
THRESH  in  chapter  6.3  could  be  adequate  in  most  cases. 


8.3.5  Test  and  results  with  procedure  LINDET 

The  objective  of  tests  was  to  identify  the  border  line 
of  lake  WALCHEN.  Border  lines  can  be  searched  either  on 
original  images  or  on  preprocessed  images.  The  two  search 
problems  differ  only  in  the  choice  of  a  different  form  of 
the  cost  function.  The  simultaneous  search  on  all  4  channels 
of  the  image  was  found  to  be  impracticable  because  of 
currently  very  long  computation  times.  Hence  the  4  channels 
have  been  compressed  according  to  equation  5.25  in  chapter 
5.5.  Parameters  were  calculated  with  the  aid  of  feature  type 
lake  (WALCHENSEE,  KOCHELSEE) .  An  image  section  with  lake 
WALCHEN  is  shown  in  figure  8.26. 


Figure  8.26;  Portion  of 
compressed  image  (Equ. 
5.25)  with  lake  WALCHEN. 
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Line  following  was  divided  into  3  segments,  which  are 
searched  independently.  The  following  cost  function  was  used. 

c(x,y)  =  255  -  |g(y)  -  g(left(y))  (  (8.4) 

where  y  is  the  pixel  on  the  right  hand  of  the  edge  and  left 
(y)  defines  its  left  hand  counterpart.  Starting  and  end 
points  of  the  three  segments  are  extremal  points  of  the 
contour  WALCHENSEE  as  designated  in  figure  8.27.  This 
segmentation  is  necessary  for  the  Euclidean  distance  used  as 
heuristic  function.  The  union  of  the  three  segments  and  the 
corresponding  development  areas  are  shown  in  figure  8.29a. 


Figure  8.2) :  Segments  of  border  line  of  WALCHENSEE, 
projection  from  map  WALCH. 


The  second  method  works  on  an  edge  image.  The  operator 
of  Roberts  was  applied  to  the  compressed  image  as  shown  in 
figure  8.28. 


Figure  8.28:  (a)  Edge  image*  of  Roberts  and  (b)  binary- 
edge  image  after  thresholding. 


Besides  the  use  of  an  Euclidean  distance  as  a  heuristic 
also  the  following  cost  function  satisfies  the  conditions 
1,2  and  3  of  chapter  6.5.2: 


cg(x)  =  (9.6859  -  0.037849  *  x)2,78  (8.5) 

The  union  of  the  three  segments  with  corresponding 
development  areas  is  shown  in  figure  8.29b. 


The  edge  representation  uses  dark  pixels  for  small  edges 
and  white  pixels  for  high  edges  except  0  edges  which  are 
transformed  to  blanks. 


Figure  8.29:  Border  line  of  lake  WALCHEN  with  development  area: 

(a)  found  in  the  compressed  original  image  and 

(b)  found  in  the  edge  image . 


Although  the  so-called  development  area  (compare, 
chapter  6.5)  in  the  second  case  (b)  is  greater,  its  completion 
time  was  less  than  in  the  first  case  (a) .  The  reason  can  be 
the  expensive  and  repetitive  computation  of  edge  values  for 
the  compressed  image,  while  for  the  edge  image  the  values  are 
already  precalculated.  This  allows  the  use  of  more  sophisticated, 
edge  operators  than  in  the  first  case. 

One  part  of  the  line  has  been  modified  at  the  left  side 
(compare  projected  object  in  figure  8.27).  The  line  follower 
skipped  the  small  branch  of  the  lake  because  it  is  not  evident 
in  the  image.  In  general  the  line  from  the  compressed  image 
is  smoother  than  the  line  from  the  edge  image  which  seems  to 
be  more  sensitive  to  little  changes  in  direction. 


8.3.6  Experiences  with  procedure  ADAPT 

As  with  procedure  THRESH,  also  ADAPT  was  applied  first 
to  object  WALCHENSEE. 


fcctfc 


The  determination  of  the  grey  level  boundary  was  done 
in  this  example  by  method  Gl  (see  chapter  8.3.4). 


A  histogram  of  the  image  in  figure  8.21  (OBI)  will 
reveal  that  pixels  are  included  that  should  be  eliminated. 

On  the  other  hand,  pixels  in  the  remainder  of  the  image  (in 
the  negative  of  OBI) ,  are  not  part  of  the  approximate  map- 
defined  feature,  but  are  adjacent  to  it  and  should  be  added. 
An  iterative  search  along  the  boundary  of  the  feature  will 
thus  permit  to  correct  the  image  (OBI) ,  using  its  histogram 
to  decide  which  pixels  to  reject  and  which  ones  to  accept. 
This 'Cinderella-algorithm* produces  the  result  shown  in 
figure  8.30a,  denoted  as  image  ABI  (adapted  binary  image). 


Figure  8.30:  Final  image  feature  after  histogram  adaption  is 

shown  in  (a) ,  whereby  the  boundary  is  marked  with 
blanks.  Figure  (b)  showes  the  differences  that  re¬ 
main  between  the  corrected  image  feature  and  the 
now  correctly  transformed  map  feature. 


Further  results  are  described  in  fi  gure  8.31  re¬ 
garding  the  second  largest  object  KOCHELSEE. 


After  these  early  experiences  from  lakes  WALCHEN  and 
KOCHEL  we  applied  ADAPT  to  seven  objects  of  the  map.  These 
objects  were  searched  on  channel  7  of  the  satellite  image 
in  two  and  in  three  iteration  steps.  Results  are  listed  in 
figure  8.32. 


Feature 

Differences 

Threshold 

boundary 

Correlation 

2.  3. 

Iteration 

2.  3. 

Iteration 

+ 

- 

WALCHENSEE 

339 

409 

(-3,10) 

(  7,30) 

95.239 

94.212 

KOCHELSEE 

181 

215 

(-3,18) 

(12,41) 

93.440 

92.176 

WA-SEE-ORT 

52 

71 

(25,46) 

(  6,31) 

82.011 

75.602 

SCHLEDORF-UF 

180 

230 

(20,44) 

(  4,37) 

71.565 

63.777 

WALD- 10 

65 

65 

(15,38) 

(  3,31) 

59.963 

57.066 

KAKOPFWALD 

35 

51 

(  6,36) 

(  0.31) 

89.668 

84.651 

KAKOPFFREI 

46 

56 

(22,48) 

(-3,23) 

56.575 

52.052 

Figure  8.32:  Results  of  method  ADAPT  on  MSS7  of  image  WALCH. 


The  threshold  boundary  was  determined  by  method  G1  (see 
chapter  8.3.4)  for  adding  pixels  (+)  working  on  a  histogram 
mask  H3 (0)  and  for  determination  of  pixels  (-)  working  on  the 
complement  of  H3 (0) .  The  resulting  grey  level  intervals  over¬ 
lap  in  some  cases  very  much,  so  that  the  same  pixel,  which 
has  been  added  in  the  first  step,  is  eliminated  in  the  second 
step.  Furthermore  no  grey  level  adaptation  was  performed  from 
iteration  to  iteration.  These  disadvantages  led  us  to  the 
procedure  formulated  in  chapter  6.4. 

Difficulties  arose  for  verification  because  with  in¬ 
creasing  iterations  the  correlation  of  the  adapted  feature 
with  the  projected  decreased  in  all  cases,  although  the 
centers  of  gravity  of  the  projection  moves  towards  the  visual¬ 
ly  perceived  collection  of  feature  pixels. With  each  iteration 
ADAPT  changes  slightly  the  shape  of  the  feature  so  that  it 
can  become  less  similar  to  the  shape  from  the  map.  This  can 
as  a  rule  be  the  case  at  the  beginning  of  ADAPT.  It  can  be 
compensated  at  higher  iterations  if  in  essence  the  image 
features  is  identical  in  shape  to  the  map  feature.  A  pre¬ 
requisite  for  reaching  an  improved  similarity  is  the  sufficient¬ 
ly  high  number  of  ADAPT-iterations . 

Figure  8.33  has  been  prepared  to  show  that  the  deterio¬ 
ration  of  the  correlation  measure  does  not  mean  a  concurrent 
motion  of  the  image  feature  into  an  erroneous  direction.  It 
compares  ADAPT  to  SHIFT  results,  listing  the  difference 
vectors  between  corresponding  centers  of  gravity.  While  cor¬ 
relation  deteriorates,  difference  vectors  do  not  change  signi¬ 
ficantly. 


Feature 


After  2  iterations 


After  3  iterations 


Difference  vector 

Length 

Difference  vector 

Length 

WALCHENSEE 

(  0.151 

9 

-0.849) 

0.863 

(  0.391 

9 

-0.713) 

0.813 

KOCHELSEE 

(-0.278 

9 

-1.004) 

1.042 

(-0.585 

9 

-1.108) 

1.253 

WA-SEE-ORT 

(-0.844 

9 

-0.699) 

1.096 

(-0.883 

9 

-0.880) 

1.247 

SCHLEDORF-UF 

(  1.518 

9 

1.824) 

2.374 

(  1.706 

9 

3.149) 

3.582 

WALD- 10 

(-1.606 

9 

-0.900) 

1.841 

(-1.695 

9 

-0.184) 

1.705 

KAKOPFWALD  - 

(-1.404 

9 

0.425) 

1.467 

(-1.327 

9 

0.461) 

1.404 

KAKOPFFREI 

(-1.725 

9 

-0.355) 

1.761 

(-1.841 

9 

-0.423) 

1.889 

Figure  8.33:  Differences  between  results  of  SHIFT  and  ADAPT 
with  2  and  3  iterations. 


For  objects  SCHLEDORF-UF  and  WALD-10  of  type  grassland  the 
automatic  grey  level  determination  found  the  wrong  intervals. 
Objects  KOCHELSEE  and  WA-SEE-ORT deteriorate  little  while  objects 
WALCHENSEE ,  WALD-10  and  KAKOPFWALD  show  decreasing  distances 
to  the  estimated  centers  of  gravity  from  procedure  SHIFT.  The 
problem  of  an  effective  verification  could  only  be  resolved 
by  allowing  a  maximum  iteration  count  much  higher  than  in  the 
present  case.  Computing  time  considerations  prevented  us  from 
more  extensive  tests. 


8.3.7  Rectification  of  image  WALCH 

Image  WALCH  was  rectified  using  13  control  points 
that  were  defined  automatically  with  method  SHIFT. 

The  13  features  served  to  compute  centers  of  gravity. 
This  led  then  to  a  set  of  control  points  for  a  rectification. 
Figures  8.34a,b  present  the  result  of  the  rectification. 


Figure  8.34:  LANDSAT- scene  from  Southern  Germany,  with  raw 
version  in  (a)  and  rectified  image  in  (b) , 
using  13  features. 


8.4  RESULTS  FROM  LANDSAT  IMAGE  GRAZ 

Whereas  image  WALCH  was  used  to  study  feature  recognition 
methods  and  their  suitability  in  connection  with  automatic 
map-image  registration,  image  GRAZ  had  the  purpose  to  test 
automation  and  performance  of  the  overall  concept. 


8.4.1  Image  GRAZ  and  map 

The  original  satellite  image  channels  4,  5,  6  and  7  of 
the  area  of  interest  chosen  for  this  test vere shown  in  figure 
5. 15 (a) -(d).  It  is  an  industrial  region  with  the  city  of  GRAZ, 
and  with  river  MUR  from  left  top  to  right  bottom.  Two  high¬ 
ways  are  marked  linear  features.  A  compressed  image  especially 
preprocessed  for  feature  type  water  (river  MUR,  see  figure 
5.18)  was  used  for  line  following.  In  the  same  way  the  four 
original  image  channels  are  compressed  for  feature  types 
buildings  (3),  grassland  (4)  and  forests  (5).  For  the  high¬ 
ways  a  ratio  image  of  channels  7  and  4  was  used.  Figure  8.35 
shows  a  printout  of  the  area  of  interst. 


Figure  8.35:  Ratio  image  (MSS7/ 

MSS 4)  of  GRAZ. 

A  plot  of  the  corresponding  map  was  shown  in  figure 
3.2  (chapter  3.2.3).  A  grid  overlay  helps  to  locate  map 
coordinates  in  that  plot.  It  should  be  stressed  that  the 
digital  map  data  bank  was  created  by  an  operator  who  did 
not  work  with  the  purpose  of  satellite  feature  recognition. 
Therefore  the  map  contains  political  boundaries  and  other 
features  that  can  be  predicted  not  to  be  identifiable  from 
a  LANDSAT  image  . 

A  list  of  areal  features  which  have  been  used  for  recogni¬ 
tion,  is  given  in  t4.g.  8.36.  This  list  is  sorted  automatically 
according  to  the  object  size  of  the  projection  in  the  image 
and  containes  further  a  point  of  the  feature  (center  of 
gravity)  in  sap  coordinates  and  following  types:  3  for 
buildings,  4  for  grassland  and  5  for  forests. 


Feature 


Number 

of 

pixel 


Type 


Center  of  gravity 
in  map  coordinates 


WALTENDORF 

3118 

3 

(225.511  ,  2046.383) 

KAISERWALD-S 

3114 

5 

(218.680  ,  2016.149) 

SCHOTTERGRUB 

1712 

4 

(219.796  ,  2031.103) 

EGGENBERG 

1429 

3 

(214.445  ,  2049.690) 

KAISERWALD-N 

1143 

5 

(213.668  ,  2021.387) 

ROSENBERG 

825 

5 

(220.635  ,  2053.956) 

KOLLISCH 

800 

5 

(231.189  ,  2012.443) 

MURFELD1 

564 

4 

(223.306  ,  2039.623) 

MESSENDORF 

538 

4 

(226.511  ,  2041.478) 

GEIDORF 

506 

3 

(221.888  ,  2051.610) 

WUNDSCHUH 

337 

5 

(224.708  ,  2019.183) 

MURFELD2 

290 

4 

(223.465  ,  2037.649) 

FORST 

258 

5 

(219.416  ,  2029.027) 

GREITH 

204 

4 

(230.174  ,  2011.791) 

PUCHWERK 

89 

3 

(226.224  ,  2036.672) 

SCHLOSSBERG 

48 

5 

(220.120  ,  2049.480) 

FLIEGERHORST 

34 

3 

(219.906  ,  2029.259) 

THALERHOF 

19 

3 

(221.295  ,  2031.384) 

Figure  8.36:  Sorted  object  list  for  image  and  map  GRAZ. 


8.4.2  Approximate  projection  from  the  map  to  the  image. 


The  projection  for  image  WALCH  was  approximated  manually 
rather  than  with'  satellite  orbit  data;  a  set  of  control 
points  was  defined  by  hand  and  a  transformation  was  calculated. 
For  image  GRAZ  the  approximate  transformation  function  was 
calculated  totally  from  informations  delivered  on  the  CCT 
magnetic  tape  header  of  the  satellite  image.  The  calculation 
was  based  on  the  satellite's  position  (47°i7')  and  attitude 
and  some  specific  properties  of  the  LANDSAT- sensor  according 
to  descriptions  of  ANUTA  (1973)  and  MALILA  et  al. (1973) .  The 
resulting  transformation  had  the  following  parameters: 


Tx<x,y)  =  8.5256  *  x  -  2.5192  *  y  +  3429.44 
Ty(x,y)  =-1.4883  *  x  -  6.2442  *  y  +  13177.08 


(8.5) 


Visual  inspection  of  some  projected  objects  verified  that  the 
accuracy  was  within  the  expected  maximum  allowance  of  about  11 
pixels . 


8.4.3  Results  with  procedure  SHIFT 

Procedure  SHIFT  was  applied  to  the  first  13  objects  of 
the  list  (figure  8.36)  sorted  by  decreasing  object  sizes. 
The  correlation  measure  was  calculated  between  the  binary 
matrix  representing  the  object  and  a  compressed  image  which 
was  preprocessed  according  to  the  feature  type. 
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The  maximum  shift  along  the  coordinate  axis  was  chosen  as  11  . 
In  accordance  with  the  papers  mentioned  above  to  calculate 
the  transformation  parameters  the  error  of  such  a  transfor¬ 
mation  can  amount  to' about  11  pixels  or  less  . 


Figure  8.37  lists  the  results.  They  are  sorted  by 
decreasing  correlation  values  so  that  the  best  fitting  of 
an  object  with  the  image  stands  on  the  top  (FORST) .  The 
resulting  relative  shift  vector  to  attain  the  maximum 
correlation  from  the  start  projection  is  used  to  calculate 
control  points  for  the  objects  with  a  correlation  greater 
than  50  X. 


Feature 

Correlation 

[X] 

Optimum  shift 
(col. , scan  lines) 

FORST 

68.667 

(-  5 

9 

BIB 

KAISERWALD-S 

66.228 

(-11 

9 

mu 

WUNDSCHUH 

64.603 

(-  8 

9 

tsi 

KAISERWALD-N 

60.159 

(-  7 

9 

mm 

KOLLISCH 

36.379 

(-  7 

9 

-i) 

MURFELD2 

34.459 

9 

+3) 

SCHOTTERGRUB 

33.834 

9 

+2) 

MESSENDORF 

30.945 

(-10 

9 

+2) 

EGGENBERG 

25.285 

(-  1 

9 

+5) 

GEIDORF 

17.403 

(-  9 

9 

+3) 

WALTENDORF 

15.957 

(-11 

9 

+2) 

ROSENBERG 

15.327 

(-  1 

9 

+2) 

MURFELD1 

9.335 

(+  3 

9 

+4) 

Figure  8.37:  Results  of  procedure  SHIFT  for  GRAZ. 


Only  four  objects  fulfil  this  verification  condition:  FORST, 
KAISERWALD-S,  WUNDSCHTJH  and  KAISERWALD-N.  These  are  all  of 
type  forest  (5).  A  subsequent  optical  inspection  of  the 
image  confirmed  the  following  experience: 

Objects  can  be  automatically  recognized  on  digital 
satellite  images  if  the  object  can  be  clearly  distinguished 
by  a  human  interpreter. 

This  experiences  is  obvious  of  course  and  is  reflected 
in  publications  on  pattern  recognition.  It  must  be  borne  in 
mind  in  the  present  purpose. 


An  important  reason  for  the  small  correlation  values 
can  be  found  in  the  map  data  bank.  WALTENDORF  for  example  is 
a  political  district,  but  was  erroneously  coded  as  a  "built-up 
area";  it  could  never  be  recognized  in  a  satellite  image. 
KOLLISCH  and  ROSENBERG  have  changed  their  shape  because  the 
map  has  not  been  revised  to  reflect  recent  urban  growth. 
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We  find  thus  that  digitization  of  current  topographic  maps 
can  lead  to  an  inappropriate  data  bank  for  automatic  feature 
recognition.  It  needs  to  be  complemented  by  human  intelligence 
to  flag  features  that  can  predictably  not  be  recognized. 

’for  comparison  reasons  the  four  recognized  objects  are 
used  to  compute  a  new  6-parameter  transformation  by  least 
squares.  Figure  8.38  lists  the  resulting  residuals  and  the 
corresponding  standard  deviations. 


Feature 

Error  vector 

Length 

Correl . 

KAISERWALD-N 

(-0.362  ,  -0.711) 

0.798 

60.159  X 

KAISERWALD-S 

(  0.466  ,  0.915) 

1.027 

66.228  X 

WUNDSCHUH 

(-0.329  ,  -0.645) 

0.724 

64.603  X 

FORST 

(  0.225  ,  0.441) 

0.495 

68.667  X 

Stand. Deviation 

(  0.712  ,  1.398) 

1.569 

64.914  X 

Figure  8.38:  Residuals  of  the  four  objects  recognized  by 
SHIFT. 


8.4.4  Results  with  procedure  THRESH 

Several  parameters  and  configurations  could  be  determined 
in  project  WALCH.  From  the  experiment  with  lake  WALCHEN,  where 
we  eliminated  pixels  of  the  pattern  after  the  threshold 
operation  according  to  their  distance  to  the  map  projection, 
we  learned  to  restrict  the  threshold  operation  to  a  starting 
mask.  This  starting  mask  has  an  equivalent  effect  of  not 
taking  into  account  pixels  far  away  from  the  projected  object. 

By  this  change  the  philosophy  of  the  grey  bound  determination 
must  be  adapted. 

The  question  arises  what  factor  p  should  be  used  (line  6  of 
procedure  THRESH,  figure  6.3a).  Since  the  resulting  pattern 
should  have  at  least  as  many  pixels  as  the  projection,  the 
factor  must  be  at  least  1.0.  For  the  upper  bound  we  empirically 
chose  1.5.  This  interval  [1.0,  1.5]  had  to  be  tested. 

Figure  8.39  lists  the  resulting  correlation  values  for 
p  values  from  1.0  to  1.5  with  step  0.05  and  for  eight  objects. 
These  objects  have  been  chosen  according  to  previous  experiences 
(forests,  optically  distinguishable  objects).  Maximum  correlation 
is  achieved  for  different  p  values,  but  correlation  values  vary 
very  little  so  that  the  objects  recognized  at  the  maximum  are 
also  recognized  for  every  p  in  the  chosen  range.  In  the  same 
way  the  corresponding  control  point  differences  remain  less  than 
1  pixel.  Therefore  the  value  p  ■  1.10  was  selected  to  be  a 
good  compromise  for  the  four  recognized  objects  (p  =  1.10 
minimizes  correlation  differences  with  least  squares) . 


p- value  1.00  1.05  1.10  1.15  1.20  1.25  1.30  1.35  1.40 
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The  centers  of  gravity  of  the  four  "recognized*  objects 
are  then  used  to  calculate  a  6-parameter  transformation. 
Results  are  shown  in  figure  8.40.  No  explanation  could  be 
found  for  the  fact  that  the  first  component  (^position  within 
a  scan  line} of  the  error  vector  is  greater  than  1  pixel  while 
the  second  component  is  always  smaller  than  1  pixel. 


Feature 

Error  vector 

Length 

Correl . 

KAISERWALD-N 

(-1.949  ,  -0.086) 

1.951 

57.155  x 

KAISERWALD-S 

(  2.508  ,  0.111) 

2.510 

78.915  X 

WUNDSCHUH 

(-1.768  ,  -0.078) 

1.770 

68.423  X 

FORST 

(  1.209  ,  0.054) 

1.211 

65.874  X 

St  and  ..Devi  at  ion 

(  3.831  ,  0.170) 

3.835 

67.592  X 

Figure  8.40:  6-parameter  transformation  with  control  points 
from  procedure  THRESH. 


8.4.5  Experiences  of  tests  with  procedures  LISU  and  LINDET 

Two  dominant  linear  features  appear  in  the  area  of 
interest : 

-  river  MUR  and 

-  two  intersecting  highways  (AUTOBAHN) 

Figure  8.41  shows  for  example  the  "mapline"  (see  figure  6.6) 
of  the  river  MUR  and  the  corresponding  starting  mask. 

Linefollowing  worked  for  river  MUR  on  the  complementary 
image  of  figure  5.18.  This  makes  sure  that  pixels  on  the 
river  get  high  values.  The  following  cost  functions  are  results 
of  a  histogram  over  the  starting  masks: 

cg(g)  =  (25.08  -  0.09796  *  g)2,21  for  MUR 

and  .  .  (8.6) 

cg(g)  =  (48.2977  -  0.1886  * g)  for  AUTOBAHN 

With  these  parameters  procedure  LISU  was  applied  to  the 
processed  images.  The  search  for  the  highways  was  split  into 
4  separate  searches  in  a  way  that  the  end  point  was  always  the 
intersection. After  the  searches  the  4  parts  of  AUTOBAHN  ("bin¬ 
line")  have  been  joined  together  for  illustration  purposes 
(figure  8.42).  The  development  area  is  shown  in  figure  8.43. 


Feature 

Length 

of 

binline 

MUR 

AUTOBAHN 

282  pixel 

480  pixel 

Surface  of 
development 
area 


Computing 

time 


278.9  sec 
71.9  sec 


Table  8.5  :  Performance  of  procedure  LISU  with  two  different 
features . 


(a)  (b) 

Figure  8.42:  Resulting  binary  matrices  of  (a)  river  MUR  and 
(b)  highways  AUTOBAHN. 


Figure  8.43:  Development  areas  of  (a)  river  MUR  and  (b)  high 
ways  AUTOBAHN. 
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The  computing  time  for  AUTOBAHN  is  the  time-sum  for  the 
four  segments.  Computing  time  in  LISU  is  not  proportional  to 
the  length  of  the  line,  it  grows  with  a  high  order  dependency 
because  set  S  of  continuation  candidate  pixels  is  growing 
with  increasing  line  length  and  must  be  inspected  completely 
for  every  step  towards  the  goal. 


For  optical  verification  the  resulting  binary  matrices 
have  been  printed  out  together  with  the  compressed  image  of 
GRAZ  (figure  8.44). 
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Computational  verification  by  subdividing  both  "mapline" 
and  "binline"  into  equal  length  segments  for  which  a  6-para¬ 
meter  transformation  is  computed  shows  that  the  position  of 
the  starting  and  end  point  has  great  importance  for  automatic 
verification  measures.  Figure  8.45  lists  the  result  (error 
vectors  and  lengths)  for  a  segment  of  AUTOBAHN  and  for  MUR 
with  5  segments  each.  Starting  and  end  points  have  been 
eliminated  according  to  experiences  from  the  synthetic  image. 


Feature 

Error  vector 

Length 

MUR 

(  1.879  ,  -0.271) 

(  0.534  ,  -0.883) 

(  0.390  ,  0.095) 

(-2.205  ,  0.852) 

1.898 

1.032 

0.401 

2.364 

AUTOBAHN 

1 

(-0.516  ,  -0.046) 
(-2.395  ,  1.062) 

(  0.472  ,  0.071) 

(  1.842  ,  -0.879) 

0.519 

2.620 

0.477 

2.041 

Stand . Deviation 

(  1.921  ,  0.836) 

2.095 

Figure  8.45:  Verification  of  linear  features  by 
segmenting . 


8.4.6  Some  results  with  procedure  ADAPT 

With  procedure  SHIFT  and  THRESH  only  4  objects  could  be 
recognized  supporting  a  starting  deviation  of  the  first 
projection  of  about  11  pixels.  This  estimated  deviation  from 
the  real  position  of  the  feature  becomes  better  if  the 
results  of  procedures  SHIFT  and  THRESH  are  used  to  calculate 
the  transformation.  Control  points  of  LINDET  have  not  been 
used  because  it  could  not  be  measured  how  good  the  match  of 
corresponding  points  has  been. 

The  resulting  transformation  has  the  following  form: 


T„(x,y)  =  8.72  *  x  -  2.08  *  y  +  2492.32 

A 

Ty(x,y)  =  -1.57  *  x  -  6.26  *  y  +  13227.4 


(8.7) 


From  the  residuals  of  the  least  squares  solution  dis¬ 
crepancies  of  about  4  to  6  pixels  have  been  found. 


Another  factor  that  leads  to  increases  of  computing  time  of 
ADAPT  was  the  surface  of  the  feature.  So  5  objects  smaller 
than  for  procedures  SHIFT  and  THRESH  were  used:  WUNDSCHUH, 
GREITH,  FLIEGERHORST,  PUCHWERK  and  THALERHOF. 

The  version  of  ADAPT  as  described  in  chapter  6.4  has 
been  applied  to  these  5  objects.  Figure  8.46  lists  the  pixel 
changes  (added  +,  eliminated  -)  for  each  iteration: 


Correl. [x] 


WUNDSCHUH 

GREITH 

FLIEGERHORST 

PUCHWERK 

57 

21 

35 

27 

23 

29 

10 

32 

4 

44 

3 

38 

72 

.031 

33  20 

19  23 

11  27 

2  31 


71.183 


62.655 


39 

0 

40 

0 

42 

13 

43 

7 

46 

7 

47 

12 

49. 

337 

□ 

0 

D 

0 

15 

0 

17 

0 

46.852 


Figure  8.46:  Adapting  5  objects  to  the  satellite  image. 


The  last  line  shows  that  the  first  3  objects  have  been 
recognized  successfully.  THALERHOF  was  too  small  to  allow  a 
correct  grey  value  determination  by  the  automatic  routine. 

The  range  was  therefore  too large  and  pixels  had  not  been 
eliminated  during  ADAPT.  A  similar  effect  accured  for  PUCH¬ 
WERK.  A  reason  why  FLIEGERHORST  could  be  verified  although 
it  is  smaller  than  PUCHWERK  was  found  in  the  bordering  forest. 
It  gives  a  good  contrast  to  the  feature  so  that  the  automatic 
grey  value  determination  delivered  a  correct  interval. 

With  three  control  points  from  ADAPT,  no  overdetermination 
was  given  for  a  6-parameter  transformation  for  ADAPT  and  there¬ 
fore  no  residuals  and  standard  deviation  could  be  calculated. 


8.4.7  Rectification  of  image  GRAZ 

With  the  aid  of  the  recognised  control  points 
procedure  RECTIFY  was  applied  to  the  four  channels  of  the 
original  satellite  image.  Discrepancies  between  control  points 
of  method  SHIFT  and  method  THRESH  have  been  eliminated  by 
using  the  result  with  the  better  verification  measure.  Figure 
8.47  can  be  regarded  as  the  final  result  of  the  automatic 
registration  process  on  GRAZ. 
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Figure  8.47:  Rectified  satellite  image  of  GRAZ 


9.  CONCLUSIONS  AND  OUTLOOK 


This  report  describes  tests  which  have  been  performed  to 
verify  a  general  concept  for  automatic  registration  of  a 
digital  map  with  a  satellite  image.  Programs  and  methods 
have  been  developed  for  this  purpose.  The  versions  of  algo¬ 
rithms  as  described  here  are  products  of  a  development  process 
and  are  not  operational.  They  only  demonstrate  the  principle 
of  our  concept  and  the  fact  that  it  works  under  certain  con¬ 
straints  . 

Future  research  should  be  concentrated  on  optimization 
of  the  general  performance  of  each  of  the  proposed  procedures 
and  the  development  of  an  appropriate  processing  system.  New 
procedures  can  then  be  included  in  the  modular  concept. 

From  the  experiences,  particularly  with  image  GRAZ,  the 
following  suggestions  can  be  derived  regarding  the  required  pro¬ 
cessing  system. 

A  processing  system  for  automated  feature  recognition: 

A  processing  system  is  required  that  implements  procedure 
REGISTER  in  more  than  one  step.  Typically  an  initial  step 
should  be  based  on  satellite  orbit  data  and  large,  distinct 
map  objects,  using  a  simple  and  fast  recognition  method.  This 
result  in  a  few  control  points  which  can  be  used,  together 
with  flight  data,  to  compute  a  better  approximation  for  the  pro¬ 
jection  function.  One  can  then  more  efficiently  employ  more 
complex  methods  and  use  smaller  and  therefore  also  more  features 

The  processing  system  must  implement  a  carefully  designed 
procedure  OBJECTSELECTION  that  incorporates  probabilities  of 
being  able  to  recognize  a  feature.  Political  districts  are  of 
course  inappropriate  for  the  purpose.  A  learning  capability 
should  be  included.  Sorting  criteria  should  be  type  of  fea¬ 
ture,  size,  contrast  with  surroundings,  predictability  of  its 
appearance  in  the  image,  scores  from  past  uses  . 

Pre-  and  postprocessing  functions : 

Preprocessing  must  be  more  adapted  to  the. used  recognition 
procedures  and  object  types.  Texture  classification  should  be 
added.  Tests  should  be  done  to  study  the  effect  of  some  scale 
reducing  preprocessing  (parsing) .  Thresholding  depends  very 
much  on  preprocessing.  Therefore  typical  features  should  be 
associated  with  characteristic  preprocessing.  We  did  test 
some  postprocessing  after  thresholding  (image  WALCH,  object 
WALCHENSEE)  by  a  distance  transform  of  the  intersection  of  the 
projection  and  threshold  masks  to  eliminate  obvious  noise. 

Post-processing  after  initial  pixel-selection  is  thus  of 
particular  interest  in  THRESH. 
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Reprogramming  of  more  operational  versions  of  recognition 
techniques : 

It  is  In  the  nature  of  research  and  development  that 
computer  programs  are  general,  flexible  and  slow.  Numerous 
ideas  exist  and  could  now  be  implemented  to  considerably 
improve  the  throughput  times  for  all  investigated  algorithms. 
The  improvements  can  be  in  the  software,  using  appropriate 
computing  strategies  (compare  brute-force  similarity  detec¬ 
tion  via  correlation  factoo  and  a  sequential  similarity 
detection  as  proposed  by  BARNEA  and  SILVERMAN,  1972) .  They 
can  be  in  hardware,  using  e.g.  pipeline  and/or  array  pro¬ 
cessing  equipment. 

To  give  just  one  example,  one  can  consider  line-fol¬ 
lowing  software  where  there  is  considerable  room  for 
improvements,  limiting  set  S.  Even  a  small  change  such  as 
eliminating  points  from  S  based  on  time  when  they  were 
accepted  in  S  can  be  of  great  value.  This  would  lead  to  an 
algorithm  similar  to  that  of  MONTANARI  (see  VANDERBRUG, 

1973)  where  backtracking  is  limited  to  the  last  three 
levels . 

Verification: 

Verification  is  very  rudimentary  at  the  moment.  Rotation 
and  scale  should  be  incorporated  in  the  verification  process 
to  yield  better  discrimination  measures  for  similarity  of 
binary  matrices.  One  possibility  can  be  to  shift  the  cen¬ 
ters  of  gravity  to  be  on  the  same  place  and  then  to  correlate 
the  two  shapes. 

If  verification  improves,  then  more  than  a  single  con¬ 
trol  point  can  be  extracted  from  a  single  object. 

In  spite  of  the  fact  that  software  is  thus  far  from 
optimized  we  can  conclude  from  the  experiences  that  ARSIM  is 
a  valid  concept  for  the  computer- supported  registration  of 
maps  and  images.  Limitations  of  the  success  are  caused  by  the 
limitation  of  experiences  and  the  map  data  bank.  We  feel 
that  the  proposed  approach  to  the  problem  is  flexible,  modu¬ 
lar  and  broad  so  that  registration  can  be  done  in  the  optimum 
way.  The  problem  is  then  to  define  this  optimum.  This,  how¬ 
ever,  is  a  question  of  experimentation  and  learning. 

Recently  it  has  become  evident,  that  the  combined  use  of 
images  and  maps  is  a  necessity  to  study  complex  environmental 
and  other  questions  •  The  ARSIM-concept  may  thus  be  of  interest 
to  a  broad  range  of  applications ,  even  where  control-point  de¬ 
finition  is  at  first  sight  not  a  problem  of  automation.  Ultima¬ 
tely  it  will  be  of  great  interest  to  g;uide  the  image  analysis 
using  map  information.  The  inherent  differences  between  maps 
and  images  can  then  be  handled  through  the  logic  that  is 
implemented  in  ARSIM. 
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